]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add LinearAlgebra::Vector derived from LinearAlgebra::ReadWriteVector and LinearAlgeb...
authorBruno Turcksin <bruno.turcksin@gmail.com>
Tue, 27 Oct 2015 14:19:04 +0000 (09:19 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Mon, 9 Nov 2015 20:49:46 +0000 (14:49 -0600)
17 files changed:
include/deal.II/lac/la_vector.h [new file with mode: 0644]
include/deal.II/lac/la_vector.templates.h [new file with mode: 0644]
include/deal.II/lac/vector_space_vector.h
source/lac/CMakeLists.txt
source/lac/la_vector.cc [new file with mode: 0644]
source/lac/la_vector.inst.in [new file with mode: 0644]
source/lac/read_write_vector.inst.in
tests/lac/la_vector_accumulation_01.cc [new file with mode: 0644]
tests/lac/la_vector_accumulation_01.output [new file with mode: 0644]
tests/lac/la_vector_add_and_dot.cc [new file with mode: 0644]
tests/lac/la_vector_add_and_dot.output [new file with mode: 0644]
tests/lac/la_vector_large_numbers.cc [new file with mode: 0644]
tests/lac/la_vector_large_numbers.output [new file with mode: 0644]
tests/lac/la_vector_norms.cc [new file with mode: 0644]
tests/lac/la_vector_norms.output [new file with mode: 0644]
tests/lac/la_vector_vector.cc [new file with mode: 0644]
tests/lac/la_vector_vector.output [new file with mode: 0644]

diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h
new file mode 100644 (file)
index 0000000..20d972f
--- /dev/null
@@ -0,0 +1,311 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C)  2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__la_vector_h
+#define dealii__la_vector_h
+
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/logstream.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/base/index_set.h>
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/vector_space_vector.h>
+
+#include <cstdio>
+#include <iostream>
+#include <cstring>
+#include <vector>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  /*! @addtogroup Vectors
+   *@{
+   */
+
+  /**
+   * Numerical vector of data. This class derives from both
+   * ::dealii::LinearAlgebra::ReadWriteVector and
+   * ::dealii::LinearAlgebra::VectorSpaceVector. As opposed to the array of the
+   * C++ standard library, this class implements an element of a vector space
+   * suitable for numerical computations.
+   *
+   * @authod Bruno Turcksin, 2015.
+   */
+  template <typename Number>
+  class Vector : public ReadWriteVector<Number>, public VectorSpaceVector<Number>
+  {
+  public:
+    /**
+     * Constructor. Create a vector of dimension zero.
+     */
+    Vector();
+
+    /**
+     * Copy constructor. Sets the dimension to that of the given vector and
+     * copies all elements.
+     */
+    Vector(const Vector<Number> &V);
+
+    /**
+     * Constructor. Set dimension to @p n and initialize all elements with
+     * zero.
+     *
+     * The constructor is made explicit to avoid accident like this:
+     * <tt>v=0;</tt>. Presumably, the user wants to set every element of the
+     * vector to zero, but instead, what happens is this call:
+     * <tt>v=Vector@<Number@>(0);</tt>, i.e. the vector is replaced by one of
+     * length zero.
+     */
+    explicit Vector(const typename ReadWriteVector<Number>::size_type n);
+
+    /**
+     * Initialize the vector with a given range of values pointed to by the
+     * iterators. This function is there in analogy to the @p std::vector class.
+     */
+    template <typename InputIterator>
+    Vector(const InputIterator first, const InputIterator last);
+
+    /**
+     * Destructor, deallocates memory.
+     */
+    ~Vector();
+
+    /**
+     * Multiply the entire vector by a fixed factor.
+     */
+    VectorSpaceVector<Number> &operator*= (const Number factor);
+
+    /**
+     * Divide the entire vector by a fixed factor.
+     */
+    VectorSpaceVector<Number> &operator/= (const Number factor);
+
+    /**
+     * Add the vector @p V to the present one.
+     */
+    VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+
+    /**
+     * Substract the vector @p V from the present one.
+     */
+    VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+
+    /**
+     * Return the scalar product of two vectors.
+     */
+    Number operator* (const VectorSpaceVector<Number> &V);
+
+    /**
+     * Simple addition of a multiple of a vector, i.e. <tt>*this += a*V</tt>.
+     */
+    void add(const Number a, const VectorSpaceVector<Number> &V);
+
+    /**
+     * Multiple addition of a multiple of a vector, i.e. <tt>*this += a*V+b*W</tt>.
+     */
+    void add(const Number a, const VectorSpaceVector<Number> &V,
+             const Number b, const VectorSpaceVector<Number> &W);
+
+    /**
+     * Scaling and simple vector addition, i.e. <tt>*this = s*(*this)+V</tt>.
+     */
+    void sadd(const Number s, const VectorSpaceVector<Number> &V);
+
+    /**
+     * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
+     * s*(*this)+a*V<tt>.
+     */
+    void sadd(const Number s, const Number a,
+              const VectorSpaceVector<Number> &V);
+
+    /**
+     * Scale each element of this vector by the corresponding element in the
+     * argument. This function is mostly meant to simulate multiplication (and
+     * immediate re-assignement) by a diagonal scaling matrix.
+     */
+    void scale(const VectorSpaceVector<Number> &scaling_factors);
+
+    /**
+     * Assignement <tt>*this = a*V</tt>.
+     */
+    void equ(const Number a, const VectorSpaceVector<Number> &V);
+
+    /**
+     * Returns the l<sub>1</sub> norm of the vector (i.e., the sum of the
+     * absolute values of all entries).
+     */
+    typename VectorSpaceVector<Number>::real_type l1_norm();
+
+    /**
+     * Returns the l<sub>2</sub> norm of the vector (i.e., the square root of
+     * the sum of the square of all entries among all processors).
+     */
+    typename VectorSpaceVector<Number>::real_type l2_norm();
+
+    /**
+     * Returns the maximum norm of the vector (i.e., the maximum absolute
+     * value among all entries and among all processors).
+     */
+    typename VectorSpaceVector<Number>::real_type linfty_norm();
+
+    /**
+     * Performs a combined operation of a vector addition and a subsequent
+     * inner product, returning the value of the inner product. In other
+     * words, the result of this function is the same as if the user called
+     * @code
+     * this->add(a, V);
+     * return_value = *this * W;
+     * @endcode
+     */
+    Number add_and_dot(const Number a,
+                       const VectorSpaceVector<Number> &V,
+                       const VectorSpaceVector<Number> &W);
+
+    /**
+     * Returns the global size of the vector, equal to the sum of the number
+     * of locally owned indices among all processors.
+     */
+    typename ReadWriteVector<Number>::size_type size() const;
+
+    /**
+     * Returns an index set that describes which elements of this vector are
+     * owned by the current processor. As a consequence, the index sets returned
+     * on different procesors if this is a distributed vector will form disjoint
+     * sets that add up to the complete index set. Obviously, if a vector is
+     * created on only one processor, then the result would satisfy
+     * @code
+     *  vec.locally_owned_elements() == complete_index_set(vec.size())
+     * @endcode
+     */
+    const dealii::IndexSet &locally_owned_elements() const;
+
+    /**
+     * Prints the vector to the output stream @p out.
+     */
+    void print(std::ostream &out,
+               const unsigned int precision=3,
+               const bool scientific=true,
+               const bool across=true) const;
+
+    /**
+     * Returns the memory consumption of this class in bytes.
+     */
+    std::size_t memory_consumption() const;
+
+  private:
+    /**
+     * Compute the L1 norm in a recursive way by dividing the vector on smaller
+     * and smaller intervals. This reduces the numerical error on large vector.
+     */
+    typename VectorSpaceVector<Number>::real_type l1_norm_recursive(unsigned int i,
+        unsigned int j);
+
+    /**
+     * Compute the L2 norm in a recursive way by dividing the vector on smaller
+     * and smaller intervals. This reduces the numerical error on large vector.
+     */
+    typename VectorSpaceVector<Number>::real_type l2_norm_recursive(unsigned int i,
+        unsigned int j);
+  };
+
+  /*@}*/
+  /*----------------------- Inline functions ----------------------------------*/
+
+  template <typename Number>
+  inline
+  Vector<Number>::Vector() {}
+
+
+
+  template <typename Number>
+  inline
+  Vector<Number>::Vector(const Vector<Number> &V)
+    :
+    ReadWriteVector<Number>(V)
+  {}
+
+
+
+  template <typename Number>
+  inline
+  Vector<Number>::Vector(const typename ReadWriteVector<Number>::size_type n)
+    :
+    ReadWriteVector<Number>(n)
+  {}
+
+
+
+  template <typename Number>
+  template <typename InputIterator>
+  inline
+  Vector<Number>::Vector(const InputIterator first, const InputIterator last)
+  {
+    this->reinit(complete_index_set(std::distance (first, last)), true);
+    std::copy(first, last, this->begin());
+  }
+
+
+
+  template <typename Number>
+  inline
+  Vector<Number>::~Vector() {}
+
+
+
+  template <typename Number>
+  inline
+  typename ReadWriteVector<Number>::size_type Vector<Number>::size() const
+  {
+    return ReadWriteVector<Number>::size();
+  }
+
+
+
+  template <typename Number>
+  inline
+  const dealii::IndexSet &Vector<Number>::locally_owned_elements() const
+  {
+    return ReadWriteVector<Number>::get_stored_elements();
+  }
+
+
+
+  template <typename Number>
+  inline
+  void Vector<Number>::print(std::ostream &out,
+                             const unsigned int precision,
+                             const bool scientific,
+                             const bool) const
+  {
+    ReadWriteVector<Number>::print(out, precision, scientific);
+  }
+
+
+
+  template <typename Number>
+  inline
+  std::size_t Vector<Number>::memory_consumption() const
+  {
+    return ReadWriteVector<Number>::memory_consumption();
+  }
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h
new file mode 100644 (file)
index 0000000..9057e54
--- /dev/null
@@ -0,0 +1,319 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C)  2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii__la_vector_templates_h
+#define dealii__la_vector_templates_h
+
+#include <deal.II/lac/la_vector.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  template <typename Number>
+  VectorSpaceVector<Number> &Vector<Number>::operator*= (const Number factor)
+  {
+    AssertIsFinite(factor);
+    for (unsigned int i=0; i<this->size(); ++i)
+      this->val[i] *= factor;
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  VectorSpaceVector<Number> &Vector<Number>::operator/= (const Number factor)
+  {
+    AssertIsFinite(factor);
+    Assert(factor!=Number(0.), ExcZero());
+    this->operator *= (Number(1.)/factor);
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  VectorSpaceVector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    Assert(down_V.size()==this->size(),
+           ExcMessage("Cannot add two vectors with different numbers of elements"));
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_V[i]);
+        this->val[i] += down_V[i];
+      }
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  VectorSpaceVector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    Assert(down_V.size()==this->size(),
+           ExcMessage("Cannot subtract two vectors with different numbers of elements"));
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_V[i]);
+        this->val[i] -= down_V[i];
+      }
+
+    return *this;
+  }
+
+
+
+  template <typename Number>
+  Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    Assert(down_V.size()==this->size(),
+           ExcMessage("Cannot compute the scalar product "
+                      "of two vectors with different numbers of elements"));
+    Number value = 0.;
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_V[i]);
+        value += this->val[i]*down_V[i];
+      }
+
+    return value;
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    AssertIsFinite(a);
+    Assert(down_V.size()==this->size(),
+           ExcMessage("Cannot add two vectors with different numbers of elements"));
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_V[i]);
+        this->val[i] += a*down_V[i];
+      }
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V,
+                           const Number b, const VectorSpaceVector<Number> &W)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    // Downcast W. If fails, throws an exception.
+    const Vector<Number> &down_W = dynamic_cast<const Vector<Number>&>(W);
+    AssertIsFinite(a);
+    Assert(down_V.size()==this->size(),
+           ExcMessage("Cannot add two vectors with different numbers of elements"));
+    AssertIsFinite(b);
+    Assert(down_W.size()==this->size(),
+           ExcMessage("Cannot add two vectors with different numbers of elements"));
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_V[i]);
+        AssertIsFinite(down_W[i]);
+        this->val[i] += a*down_V[i]+b*down_W[i];
+      }
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::sadd(const Number s, const VectorSpaceVector<Number> &V)
+  {
+    *this *= s;
+    *this += V;
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::sadd(const Number s, const Number a,
+                            const VectorSpaceVector<Number> &V)
+  {
+    *this *= s;
+    // Downcast V. It fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    Vector<Number> tmp(down_V);
+    tmp *= a;
+    *this += tmp;
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
+  {
+    // Downcast scaling_factors. If fails, throws an exception.
+    const Vector<Number> &down_scaling_factors =
+      dynamic_cast<const Vector<Number>&>(scaling_factors);
+    Assert(down_scaling_factors.size()==this->size(),
+           ExcMessage("Cannot add two vectors with different numbers of elements"));
+    for (unsigned int i=0; i<this->size(); ++i)
+      {
+        AssertIsFinite(down_scaling_factors[i]);
+        this->val[i] *= down_scaling_factors[i];
+      }
+  }
+
+
+
+  template <typename Number>
+  void Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
+  {
+    // Downcast V. If fails, throws an exception.
+    const Vector<Number> &down_V = dynamic_cast<const Vector<Number>&>(V);
+    *this = down_V;
+    *this *= a;
+  }
+
+
+
+  template <typename Number>
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm()
+  {
+    Assert (this->size(), ExcEmptyObject());
+
+    typename ReadWriteVector<Number>::real_type norm = 0.;
+    norm = l1_norm_recursive(0,this->size()-1);
+
+    return norm;
+  }
+
+
+
+
+  template <typename Number>
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm()
+  {
+    Assert (this->size(), ExcEmptyObject());
+
+    // if l2_norm()^2 is finite and non-zero, the answer is computed as
+    // std::sqrt(norm_sqr()). If norm_sqrt() is infinite or zero, the l2 norm
+    // might still be finite. In that case, recompute ig (this is a rare case,
+    // so working on the vector twice is uncritical and paid off by the extended
+    // precision) using the BLAS approach with a weight, see e.g. dnrm2.f.
+    typedef typename ReadWriteVector<Number>::real_type real_type;
+    real_type norm_square = 0.;
+    norm_square = l2_norm_recursive(0,this->size()-1);
+    if (numbers::is_finite(norm_square) &&
+        norm_square>=std::numeric_limits<real_type>::min())
+      return std::sqrt(norm_square);
+    else
+      {
+        real_type scale = 0.;
+        real_type sum = 1.;
+        for (unsigned int i=0; i<this->size(); ++i)
+          {
+            if (this->val[i] != Number())
+              {
+                const real_type abs_x =
+                  numbers::NumberTraits<Number>::abs(this->val[i]);
+                if (scale < abs_x)
+                  {
+                    sum = 1. + sum * (scale/abs_x) * (scale/abs_x);
+                    scale = abs_x;
+                  }
+                else
+                  sum += (abs_x/scale) * (abs_x/scale);
+              }
+          }
+        AssertIsFinite(scale*std::sqrt(sum));
+        return scale*std::sqrt(sum);
+      }
+  }
+
+
+
+  template <typename Number>
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm()
+  {
+    typename ReadWriteVector<Number>::real_type norm = 0.;
+    for (unsigned int i=0; i<this->size(); ++i)
+      if (std::abs(this->val[i])>norm)
+        norm = std::abs(this->val[i]);
+
+    return norm;
+  }
+
+
+
+  template <typename Number>
+  Number Vector<Number>::add_and_dot(const Number a,
+                                     const VectorSpaceVector<Number> &V,
+                                     const VectorSpaceVector<Number> &W)
+  {
+    this->add(a,V);
+
+    return *this*W;
+  }
+
+
+
+  template <typename Number>
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(unsigned int i,
+      unsigned int j)
+  {
+    Assert(j>=i, ExcInternalError());
+    typename ReadWriteVector<Number>::real_type norm = 0.;
+
+    if ((j-i)!=0)
+      {
+        norm += l1_norm_recursive(i,(i+j)/2);
+        norm += l1_norm_recursive((i+j)/2+1,j);
+      }
+    else
+      norm += std::abs(this->val[i]);
+
+    return norm;
+  }
+
+
+
+  template <typename Number>
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_recursive(unsigned int i,
+      unsigned int j)
+  {
+    Assert(j>=i, ExcInternalError());
+
+    typename ReadWriteVector<Number>::real_type norm = 0.;
+    if ((j-i)!=0)
+      {
+        norm += l2_norm_recursive(i,(i+j)/2);
+        norm += l2_norm_recursive((i+j)/2+1,j);
+      }
+    else
+      norm += std::pow(std::abs(this->val[i]),2);
+
+    return norm;
+  }
+}
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
index 47184f60e0554bf44ea94ad6e94246304e6ecb7a..60e47ee77164a60fd45645841b4108b65c998694 100644 (file)
 #include <deal.II/base/mpi.h>
 #include <deal.II/base/numbers.h>
 
+
 DEAL_II_NAMESPACE_OPEN
 
+class IndexSet;
+
 namespace LinearAlgebra
 {
   template <typename Number>
   class ReadWriteVector;
 
-  class IndexSet;
 
   /*! @addtogroup Vectors
    *@{
@@ -49,6 +51,7 @@ namespace LinearAlgebra
     typedef types::global_dof_index                           size_type;
     typedef typename numbers::NumberTraits<Number>::real_type real_type;
 
+
     /**
      * Multiply the entire vector by a fixed factor.
      */
@@ -91,7 +94,7 @@ namespace LinearAlgebra
     virtual void sadd(const Number s, const VectorSpaceVector<Number> &V) = 0;
 
     /**
-     * Scaling and simple additiion of a multiple of a vector, i.e. <tt>*this =
+     * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
      * s*(*this)+a*V</tt>.
      */
     virtual void sadd(const Number s, const Number a,
@@ -154,7 +157,7 @@ namespace LinearAlgebra
     virtual size_type size() const = 0;
 
     /**
-     * Return an index set that describes which elements of this vector are
+     * Returns an index set that describes which elements of this vector are
      * owned by the current processor. As a consequence, the index sets returned
      * on different procesors if this is a distributed vector will form disjoint
      * sets that add up to the complete index set. Obviously, if a vector is
@@ -163,7 +166,7 @@ namespace LinearAlgebra
      *  vec.locally_owned_elements() == complete_index_set(vec.size())
      * @endcode
      */
-    virtual IndexSet locally_owned_elements() const = 0;
+    virtual const dealii::IndexSet &locally_owned_elements() const = 0;
 
     /**
      * Prints the vector to the output stream @p out.
@@ -178,7 +181,6 @@ namespace LinearAlgebra
      */
     virtual std::size_t memory_consumption() const = 0;
   };
-
   /*@}*/
 }
 
index b4de9e64e341857495b79a0ce088ae22db85e6d4..648edcfc282bb2ea1b1c109700974a5b0c61bb33 100644 (file)
@@ -27,6 +27,7 @@ SET(_src
   constraint_matrix.cc
   full_matrix.cc
   lapack_full_matrix.cc
+  la_vector.cc
   matrix_lib.cc
   matrix_out.cc
   parallel_vector.cc
@@ -60,6 +61,7 @@ SET(_inst
   constraint_matrix.inst.in
   full_matrix.inst.in
   lapack_full_matrix.inst.in
+  la_vector.inst.in
   parallel_vector.inst.in
   precondition_block.inst.in
   relaxation_block.inst.in
diff --git a/source/lac/la_vector.cc b/source/lac/la_vector.cc
new file mode 100644 (file)
index 0000000..ea6500a
--- /dev/null
@@ -0,0 +1,25 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/la_vector.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+#include "la_vector.inst"
+}
+
+DEAL_II_NAMESPACE_CLOSE
diff --git a/source/lac/la_vector.inst.in b/source/lac/la_vector.inst.in
new file mode 100644 (file)
index 0000000..5402820
--- /dev/null
@@ -0,0 +1,29 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1999 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+for (SCALAR : REAL_SCALARS)
+  {
+    template class Vector<SCALAR>;
+  }
+
+
+
+for (SCALAR : COMPLEX_SCALARS)
+  {
+    template class Vector<SCALAR>;
+  }
+
index eceace6fd43a5181a0163e330d0e74626a63b863..dd83ecdbca6461bb5c9bed637cd300d1e3e3ec8e 100644 (file)
@@ -32,3 +32,21 @@ for (S1, S2 : REAL_SCALARS)
   \}
 }
 
+
+for (SCALAR : COMPLEX_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    template class ReadWriteVector<SCALAR>;
+  \}
+}
+
+for (S1, S2 : COMPLEX_SCALARS)
+{
+  namespace LinearAlgebra
+  \{
+    template void ReadWriteVector<S1>::reinit<S2> (const ReadWriteVector<S2>&,
+                                          const bool);
+  \}
+}
+
diff --git a/tests/lac/la_vector_accumulation_01.cc b/tests/lac/la_vector_accumulation_01.cc
new file mode 100644 (file)
index 0000000..a4c00f8
--- /dev/null
@@ -0,0 +1,67 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check that the l2 norm is exactly the same for many runs on random vector
+// size with random vector entries. this is to ensure that the result is
+// reproducible also when parallel evaluation is done
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check_norms ()
+{
+  for (unsigned int test=0; test<20; ++test)
+    {
+      const unsigned int size = Testing::rand() % 100000;
+      LinearAlgebra::Vector<number> vec (size);
+      for (unsigned int i=0; i<size; ++i)
+        vec(i) = static_cast<number>(Testing::rand())/static_cast<number>(RAND_MAX);
+      const typename LinearAlgebra::ReadWriteVector<number>::real_type norm = vec.l2_norm();
+      for (unsigned int i=0; i<30; ++i)
+        AssertThrow (vec.l2_norm() == norm, ExcInternalError());
+
+      LinearAlgebra::Vector<number> vec2 (vec);
+      for (unsigned int i=0; i<10; ++i)
+        AssertThrow (vec2.l2_norm() == norm, ExcInternalError());
+    }
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check_norms<float>();
+  check_norms<double>();
+  check_norms<long double>();
+  check_norms<std::complex<double> >();
+  deallog << "OK" << std::endl;
+}
+
+
diff --git a/tests/lac/la_vector_accumulation_01.output b/tests/lac/la_vector_accumulation_01.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/lac/la_vector_add_and_dot.cc b/tests/lac/la_vector_add_and_dot.cc
new file mode 100644 (file)
index 0000000..4441163
--- /dev/null
@@ -0,0 +1,83 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2014 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// check that LinearAlgebra::Vector::add_and_dot works correctly
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check ()
+{
+  for (unsigned int test=0; test<5; ++test)
+    {
+      const unsigned int size = 17 + test*1101;
+      LinearAlgebra::Vector<number> v1 (size), v2(size), v3(size), check(size);
+      for (unsigned int i=0; i<size; ++i)
+        {
+          v1(i) = 0.1 + 0.005 * i;
+          v2(i) = -5.2 + 0.18 * i;
+          v3(i) = 3.14159 + 2.7183/(1.+i);
+        }
+      check = v1;
+      const number factor = 0.01432;
+
+      v1.add(factor, v2);
+      const number prod = v1 * v3;
+      const number prod_check = check.add_and_dot(factor, v2, v3);
+      if (test == 0 && types_are_equal<number,double>::value)
+        {
+          deallog << "Vector add reference:   ";
+          for (unsigned int i=0; i<size; ++i)
+            deallog << v1[i] << " ";
+          deallog << std::endl;
+          deallog << "Vector check reference: ";
+          for (unsigned int i=0; i<size; ++i)
+            deallog << check[i] << " ";
+          deallog << std::endl;
+        }
+
+      deallog << "Add and dot should be " << prod/static_cast<number>(size)
+              << ", is " << prod_check/static_cast<number>(size)
+              << std::endl;
+    }
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check<float>();
+  check<double>();
+  check<long double>();
+  check<std::complex<double> >();
+  deallog << "OK" << std::endl;
+}
+
+
diff --git a/tests/lac/la_vector_add_and_dot.output b/tests/lac/la_vector_add_and_dot.output
new file mode 100644 (file)
index 0000000..329e9f8
--- /dev/null
@@ -0,0 +1,24 @@
+
+DEAL::Add and dot should be 0.30, is 0.30
+DEAL::Add and dot should be 13.40, is 13.40
+DEAL::Add and dot should be 26.50, is 26.50
+DEAL::Add and dot should be 39.61, is 39.61
+DEAL::Add and dot should be 52.71, is 52.71
+DEAL::Vector add reference:   0.03   0.03   0.04   0.05   0.06   0.06   0.07   0.08   0.09   0.09   0.10   0.11   0.12   0.12   0.13   0.14   0.15   
+DEAL::Vector check reference: 0.03   0.03   0.04   0.05   0.06   0.06   0.07   0.08   0.09   0.09   0.10   0.11   0.12   0.12   0.13   0.14   0.15   
+DEAL::Add and dot should be 0.30, is 0.30
+DEAL::Add and dot should be 13.40, is 13.40
+DEAL::Add and dot should be 26.50, is 26.50
+DEAL::Add and dot should be 39.61, is 39.61
+DEAL::Add and dot should be 52.71, is 52.71
+DEAL::Add and dot should be 0.30, is 0.30
+DEAL::Add and dot should be 13.40, is 13.40
+DEAL::Add and dot should be 26.50, is 26.50
+DEAL::Add and dot should be 39.61, is 39.61
+DEAL::Add and dot should be 52.71, is 52.71
+DEAL::Add and dot should be (0.30,0.00), is (0.30,0.00)
+DEAL::Add and dot should be (13.40,0.00), is (13.40,0.00)
+DEAL::Add and dot should be (26.50,0.00), is (26.50,0.00)
+DEAL::Add and dot should be (39.61,0.00), is (39.61,0.00)
+DEAL::Add and dot should be (52.71,0.00), is (52.71,0.00)
+DEAL::OK
diff --git a/tests/lac/la_vector_large_numbers.cc b/tests/lac/la_vector_large_numbers.cc
new file mode 100644 (file)
index 0000000..884e053
--- /dev/null
@@ -0,0 +1,81 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+#include <limits>
+
+
+
+void check_large_numbers()
+{
+  LinearAlgebra::Vector<float> v(10);
+  v(0) = 1e13;
+  v(1) = 1e19;
+  v(2) = 1e21;
+  v(4) = 1.;
+  v(7) = 1e20;
+
+  const double correct = std::sqrt(1e26 + 1e38 + 1e42 + 1e40 + 1.);
+  AssertThrow (std::abs(v.l2_norm() - correct) < 1e-6*correct, ExcInternalError());
+
+  for (unsigned int i=0; i<v.size(); ++i)
+    v[i] = (float) 0.;
+  v(5) = 1e-30;
+  v(9) = 1e-32;
+  const double correct2 = std::sqrt(1e-60 + 1e-64);
+  AssertThrow (std::abs(v.l2_norm() - correct2) < 1e-6*correct2, ExcInternalError());
+
+  LinearAlgebra::Vector<double> w(7);
+  w(0) = 1e232;
+  w(6) = 1e231;
+  w(3) = 1e200;
+  const double correct3 = std::sqrt(100. + 1.) * 1e231;
+  AssertThrow (std::abs(w.l2_norm() - correct3) < 1e-13*correct3, ExcInternalError());
+
+  for (unsigned int i=0; i<w.size(); ++i)
+    w[i] = (float) 0.;
+  w(1) = 1e-302;
+  w(2) = 1e-303;
+  w(3) = 2e-303;
+  w(4) = 3e-303;
+  w(5) = -1e-303;
+  const double correct4 = std::sqrt(100. + 1. + 4. + 9 + 1.) * 1e-303;
+  AssertThrow (std::abs(w.l2_norm() - correct4) <= 1e-13*correct4, ExcInternalError());
+
+  deallog << "OK" << std::endl;
+}
+
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check_large_numbers();
+}
+
+
diff --git a/tests/lac/la_vector_large_numbers.output b/tests/lac/la_vector_large_numbers.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/lac/la_vector_norms.cc b/tests/lac/la_vector_norms.cc
new file mode 100644 (file)
index 0000000..7f2faad
--- /dev/null
@@ -0,0 +1,119 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+// checks that the l1, l2, lp norm is computed correctly for deal.II vectors
+// (check the summation algorithm), including an accuracy test (should not
+// lose more than 1 decimal also for 200000 vector entries)
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check_norms ()
+{
+  const number acc = 1e1*std::numeric_limits<number>::epsilon();
+  unsigned int skip = 73;
+  for (unsigned int size=1; size<200000; size+=skip)
+    {
+      // test correctness
+      if (size > 10000)
+        skip += 17;
+      LinearAlgebra::Vector<number> vec(size);
+      for (unsigned int i=0; i<size; ++i)
+        vec(i) = i+1;
+
+      const number l1_norm = vec.l1_norm();
+      AssertThrow (std::abs(l1_norm-0.5*size*(size+1)) < acc*0.5*size*(size+1),
+                   ExcInternalError());
+
+      // test accuracy of summation
+      const long double value = 3.14159265358979323846;
+      for (unsigned int i=0; i<size; ++i)
+        vec[i] = (number)value;
+      const number l1_norma = vec.l1_norm();
+      std::cout<<l1_norma<<std::endl;
+      AssertThrow (std::abs(l1_norma-value*size) < acc*size*value,
+                   ExcInternalError());
+      const number l2_norma = vec.l2_norm();
+      AssertThrow (std::abs(l2_norma-value*std::sqrt((number)size)) < acc*std::sqrt(size)*value,
+                   ExcInternalError());
+    }
+}
+
+
+template <typename number>
+void check_complex_norms ()
+{
+  const number acc = 1e1*std::numeric_limits<number>::epsilon();
+  unsigned int skip = 73;
+  for (unsigned int size=1; size<100000; size+=skip)
+    {
+      // test correctness
+      if (size > 10000)
+        skip += 17;
+      LinearAlgebra::Vector<std::complex<number> > vec(size);
+      long double sum = 0.;
+      for (unsigned int i=0; i<size; ++i)
+        {
+          vec(i) = std::complex<number>(i+1,i+2);
+          sum += std::sqrt((long double)(i+1)*(i+1) + (long double)(i+2)*(i+2));
+        }
+
+      const number l1_norm = vec.l1_norm();
+      AssertThrow (std::abs(l1_norm-sum) < acc*sum,
+                   ExcInternalError());
+
+      // test accuracy of summation
+      const std::complex<long double> value (3.14159265358979323846, 0.1);
+      for (unsigned int i=0; i<size; ++i)
+        vec[i] = std::complex<number>(value);
+      const number l1_norma = vec.l1_norm();
+      AssertThrow (std::abs(l1_norma-std::abs(value)*size) <
+                   acc*size*std::abs(value),
+                   ExcInternalError());
+      const number l2_norma = vec.l2_norm();
+      AssertThrow (std::abs(l2_norma-std::abs(value)*std::sqrt((number)size)) <
+                   acc*std::sqrt((number)size)*std::abs(value),
+                   ExcInternalError());
+    }
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  check_norms<float>();
+  check_norms<double>();
+  check_norms<long double>();
+  check_complex_norms<double>();
+  check_complex_norms<float>();
+  deallog << "OK" << std::endl;
+}
+
+
diff --git a/tests/lac/la_vector_norms.output b/tests/lac/la_vector_norms.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/lac/la_vector_vector.cc b/tests/lac/la_vector_vector.cc
new file mode 100644 (file)
index 0000000..1dbc2c4
--- /dev/null
@@ -0,0 +1,159 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 1998 - 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/la_vector.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+const unsigned int N=10;
+unsigned int check_point = 0;
+
+
+
+
+template <typename number>
+void print (const LinearAlgebra::Vector<number> &v)
+{
+  for (unsigned int i=0; i<v.size(); ++i)
+    deallog << v(i) << '\t';
+  deallog << std::endl;
+}
+
+
+
+template <typename number1, typename number2>
+void check_vectors (LinearAlgebra::Vector<number1> &d1, LinearAlgebra::Vector<number2> &d2)
+{
+  deallog << "Fill & Swap" << std::endl;
+  LinearAlgebra::Vector<number1> d3(d1.size());
+  print (d3);
+
+  for (unsigned int i=0; i<N; ++i)
+    {
+      d1(i) = 2. * i;
+      d2(i) = .5 * d1(i)*d1(i);
+      d3(i) = 2. - .5 * i;
+    };
+
+  print (d1);
+  print (d2);
+  print (d3);
+
+  swap (d1, d2);
+  print (d1);
+
+  d1 = d2;
+  print (d1);
+
+  for (unsigned int i=0; i<N; ++i)
+    d1[i] = 2.5;
+  print (d1);
+
+  // initialize with iterators
+  number1 array[] = { 0.0, 1.1, 2.2, 3.3 };
+  LinearAlgebra::Vector<number1> d4 (&array[0], &array[4]);
+  print (d4);
+
+  deallog << "Extract number" << std::endl;
+  // Each line should contain two equal numbers
+  double sum = 0.;
+  for (unsigned int i=0; i<N; ++i)
+    sum += 4.*i-i*i;
+  deallog << d3 *d2 << '\t' << sum << std::endl;
+
+  sum = 0.;
+  for (unsigned int i=0; i<N; ++i)
+    sum += 4.*i*i;
+  sum = std::sqrt(sum);
+  deallog << d2.l2_norm() << '\t' << sum << std::endl;
+
+  sum = 0.;
+  for (unsigned int i=0; i<N; ++i)
+    sum += std::fabs(2.-.5*i);
+  deallog << d3.l1_norm() << '\t' << sum << std::endl;
+
+  sum = 0.;
+  for (unsigned int i=0; i<N; ++i)
+    {
+      double t = std::fabs(2.-.5*i);
+      if (t>sum) sum = t;
+    }
+  deallog << d3.linfty_norm() << '\t' << sum << std::endl;
+
+  deallog << "add & sub" << std::endl;
+
+  d1 += d2;
+  print (d1);
+
+  d2 -= d1;
+  print (d2);
+
+  d1.add (2, d3);
+  print (d1);
+
+  d1.add (2., d2, .5, d3);
+  print (d1);
+
+  deallog << "sadd & scale" << std::endl;
+
+  d2.sadd (2., d1);
+  print (d2);
+
+  d2.sadd (2., .5, d1);
+  print (d2);
+
+  d1 *= 4.;
+  print (d1);
+
+  deallog << "equ" << std::endl;
+
+  d2.equ (.25, d1);
+  print (d2);
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::fixed;
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1.e-10);
+
+  LinearAlgebra::Vector<double>      d1(N), d2(N);
+  LinearAlgebra::Vector<float>       f1(N), f2(N);
+  LinearAlgebra::Vector<long double> l1(N), l2(N);
+
+  // cross-tests with double/float
+  // vectors at the same time are not
+  // supported at present,
+  // unfortunately, as many functions
+  // don't accept other data types as
+  // arguments
+  check_vectors (d1, d2);
+  check_vectors (f1, f2);
+  check_vectors (l1, l2);
+}
+
+
diff --git a/tests/lac/la_vector_vector.output b/tests/lac/la_vector_vector.output
new file mode 100644 (file)
index 0000000..46bd087
--- /dev/null
@@ -0,0 +1,76 @@
+
+DEAL::Fill & Swap
+DEAL::0        0       0       0       0       0       0       0       0       0       
+DEAL::0        2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::0        2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::2.00     1.50    1.00    0.50    0       -0.50   -1.00   -1.50   -2.00   -2.50   
+DEAL::0        2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::0        2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::2.50     2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    
+DEAL::0        1.10    2.20    3.30    
+DEAL::Extract number
+DEAL::-105.00  -105.00
+DEAL::33.76    33.76
+DEAL::12.50    12.50
+DEAL::2.50     2.50
+DEAL::add & sub
+DEAL::2.50     4.50    6.50    8.50    10.50   12.50   14.50   16.50   18.50   20.50   
+DEAL::-2.50    -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   
+DEAL::6.50     7.50    8.50    9.50    10.50   11.50   12.50   13.50   14.50   15.50   
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    
+DEAL::sadd & scale
+DEAL::-2.50    -1.75   -1.00   -0.25   0.50    1.25    2.00    2.75    3.50    4.25    
+DEAL::-3.75    -1.88   0       1.88    3.75    5.62    7.50    9.38    11.25   13.12   
+DEAL::10.00    13.00   16.00   19.00   22.00   25.00   28.00   31.00   34.00   37.00   
+DEAL::equ
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    
+DEAL::Fill & Swap
+DEAL::0.00     0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    
+DEAL::0.00     2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::0.00     2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::2.00     1.50    1.00    0.50    0.00    -0.50   -1.00   -1.50   -2.00   -2.50   
+DEAL::0.00     2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::0.00     2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::2.50     2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    
+DEAL::0.00     1.10    2.20    3.30    
+DEAL::Extract number
+DEAL::-105.00  -105.00
+DEAL::33.76    33.76
+DEAL::12.50    12.50
+DEAL::2.50     2.50
+DEAL::add & sub
+DEAL::2.50     4.50    6.50    8.50    10.50   12.50   14.50   16.50   18.50   20.50   
+DEAL::-2.50    -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   
+DEAL::6.50     7.50    8.50    9.50    10.50   11.50   12.50   13.50   14.50   15.50   
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    
+DEAL::sadd & scale
+DEAL::-2.50    -1.75   -1.00   -0.25   0.50    1.25    2.00    2.75    3.50    4.25    
+DEAL::-3.75    -1.88   0.00    1.88    3.75    5.62    7.50    9.38    11.25   13.12   
+DEAL::10.00    13.00   16.00   19.00   22.00   25.00   28.00   31.00   34.00   37.00   
+DEAL::equ
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    
+DEAL::Fill & Swap
+DEAL::0.00     0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    0.00    
+DEAL::0.00     2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::0.00     2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::2.00     1.50    1.00    0.50    0.00    -0.50   -1.00   -1.50   -2.00   -2.50   
+DEAL::0.00     2.00    8.00    18.00   32.00   50.00   72.00   98.00   128.00  162.00  
+DEAL::0.00     2.00    4.00    6.00    8.00    10.00   12.00   14.00   16.00   18.00   
+DEAL::2.50     2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    2.50    
+DEAL::0.00     1.10    2.20    3.30    
+DEAL::Extract number
+DEAL::-105.00  -105.00
+DEAL::33.76    33.76
+DEAL::12.50    12.50
+DEAL::2.50     2.50
+DEAL::add & sub
+DEAL::2.50     4.50    6.50    8.50    10.50   12.50   14.50   16.50   18.50   20.50   
+DEAL::-2.50    -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   -2.50   
+DEAL::6.50     7.50    8.50    9.50    10.50   11.50   12.50   13.50   14.50   15.50   
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    
+DEAL::sadd & scale
+DEAL::-2.50    -1.75   -1.00   -0.25   0.50    1.25    2.00    2.75    3.50    4.25    
+DEAL::-3.75    -1.88   0.00    1.88    3.75    5.62    7.50    9.38    11.25   13.12   
+DEAL::10.00    13.00   16.00   19.00   22.00   25.00   28.00   31.00   34.00   37.00   
+DEAL::equ
+DEAL::2.50     3.25    4.00    4.75    5.50    6.25    7.00    7.75    8.50    9.25    

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.