]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Combined add and inner product for vectors: add_and_dot
authorMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Sun, 26 Oct 2014 08:25:43 +0000 (09:25 +0100)
committerMartin Kronbichler <kronbichler@lnm.mw.tum.de>
Fri, 31 Oct 2014 18:43:44 +0000 (12:43 -0600)
- Implemented a new method add_and_dot for all vector classes (deal.II's vectors,
  block vectors, PETSc, Trilinos vectors)
- Specialized function for deal.II vector that does add_and_dot in one sweep over
  the data by a new AddAndDot struct that is passed to the accumulate function. This
  also extends the other structs for inner products and norms by one (unused) argument.
- New tests for this feature.
- Added new do_vectorized operation to all accumulation operations that uses the
  operations of VectorizedArray<double/float> for the sums. This helps to reduce the
  operations and thus increase performance when memory transfer is not the limit.
  This functionality is called in the regular part of the accumulate function which
  now sits in its own method 'accumulate_regular' (length divisible by 32). This
  functionality can unfortunately not be realized by the OPENMP_SIMD pragma recently
  introduced because that leads to non-reproducible results when the vector memory
  starts at different memory locations. VectorizedArray is only supported on some
  platforms and thus this functionality may not be available everywhere. However, good
  code will be generated unconditionally (except on old machines from around 2004-2008
  where unaligned reads also from aligned memory are slow, but those are not relevant
  any more).
- To implemented vectorized pow operations, I also added a std::pow function for
  VectorizedArray.
- The memory for parallel::distributed::Vector is now allocated to 64 byte boundaries
  just as is the memory of Vector<Number>.

29 files changed:
doc/news/changes.h
include/deal.II/base/vectorization.h
include/deal.II/lac/block_vector_base.h
include/deal.II/lac/parallel_block_vector.h
include/deal.II/lac/parallel_vector.h
include/deal.II/lac/parallel_vector.templates.h
include/deal.II/lac/petsc_vector_base.h
include/deal.II/lac/solver_bicgstab.h
include/deal.II/lac/solver_cg.h
include/deal.II/lac/solver_gmres.h
include/deal.II/lac/trilinos_vector_base.h
include/deal.II/lac/vector.h
include/deal.II/lac/vector.templates.h
source/lac/petsc_vector_base.cc
tests/lac/gmres_reorthogonalize_02.cc
tests/lac/gmres_reorthogonalize_02.compiler=Intel=no.release.output [deleted file]
tests/lac/gmres_reorthogonalize_02.compiler=Intel=yes.release.output [deleted file]
tests/lac/gmres_reorthogonalize_02.debug.output [deleted file]
tests/lac/gmres_reorthogonalize_02.output [new file with mode: 0644]
tests/lac/gmres_reorthogonalize_04.cc
tests/lac/gmres_reorthogonalize_04.compiler=Intel=no.release.output [deleted file]
tests/lac/gmres_reorthogonalize_04.debug.output [deleted file]
tests/lac/gmres_reorthogonalize_04.output [moved from tests/lac/gmres_reorthogonalize_04.compiler=Intel=yes.release.output with 70% similarity]
tests/lac/vector_accumulation_01.cc [moved from tests/lac/vector_accumulation.cc with 100% similarity]
tests/lac/vector_accumulation_01.output [moved from tests/lac/vector_accumulation.output with 100% similarity]
tests/lac/vector_accumulation_02.cc [new file with mode: 0644]
tests/lac/vector_accumulation_02.output [new file with mode: 0644]
tests/lac/vector_add_and_dot.cc [new file with mode: 0644]
tests/lac/vector_add_and_dot.output [new file with mode: 0644]

index ce74a29e2bad79453a4d3e3f8cae3b1ed14e44d2..eb0567c0d33624ca0c9bb730d4cfb87f00f95f58 100644 (file)
@@ -307,6 +307,21 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li> New: The vector classes in deal.II (including Trilinos and PETSc
+  wrappers) now have a new method x.add_and_dot(factor,v,w) which performs
+  x.add(factor,v) and subsequent inner product of x with another vector
+  w. This operation occurs in some iterative solvers; by a combined operation,
+  reduced memory transfer and thus higher performance are enabled.
+  <br>
+  (Martin Kronbichler, 2014/10/27)
+  </li>
+
+  <li> Improved: Inner products and norms on deal.II's own vector classes now
+  use vectorization through VectorizedArray if available.
+  <br>
+  (Martin Kronbichler, 2014/10/27)
+  </li>
+
   <li> Changed: PETSc and Trilinos vectors with ghost entries can now be reset to zero
   using = 0.0;
   <br>
index d0993978b9b25cc0b8ad42af47c80f7ba03ed44c..7b903e19dbd4e81b665491541bb5ad791050c887 100644 (file)
@@ -2016,6 +2016,29 @@ namespace std
 
 
 
+  /**
+   * Raises the given number @p x to the power @p p for a vectorized data
+   * field. The result is returned as vectorized array in the form
+   * <tt>{pow(x[0],p), pow(x[1],p), ..., pow(x[n_array_elements-1],p)}</tt>.
+   *
+   * @relates VectorizedArray
+   */
+  template <typename Number>
+  inline
+  ::dealii::VectorizedArray<Number>
+  pow (const ::dealii::VectorizedArray<Number> &x,
+       const Number p)
+  {
+    Number values[::dealii::VectorizedArray<Number>::n_array_elements];
+    for (unsigned int i=0; i<dealii::VectorizedArray<Number>::n_array_elements; ++i)
+      values[i] = std::pow(x[i], p);
+    ::dealii::VectorizedArray<Number> out;
+    out.load(&values[0]);
+    return out;
+  }
+
+
+
   /**
    * Computes the absolute value (modulus) of a vectorized data field. The
    * result is returned as vectorized array in the form <tt>{abs(x[0]),
index 2a8f34580878eb7b96ac2152475efd407f516d2c..543b4a495ee6e707bb21486360c0b5dbc72fc005 100644 (file)
@@ -1024,6 +1024,27 @@ public:
    */
   real_type linfty_norm () const;
 
+  /**
+   * 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
+   *
+   * The reason this function exists is that this operation involves less
+   * memory transfer than calling the two functions separately on deal.II's
+   * vector classes (Vector<Number> and parallel::distributed::Vector<double>).
+   * This method only needs to load three vectors, @p this, @p V, @p W,
+   * whereas calling separate methods means to load the calling vector @p this
+   * twice. Since most vector operations are memory transfer limited, this
+   * reduces the time by 25\% (or 50\% if @p W equals @p this).
+   */
+  value_type add_and_dot (const value_type       a,
+                          const BlockVectorBase &V,
+                          const BlockVectorBase &W);
+
   /**
    * Returns true if the given global index is
    * in the local range of this processor.
@@ -2091,6 +2112,25 @@ BlockVectorBase<VectorType>::linfty_norm () const
 
 
 
+template <class VectorType>
+typename BlockVectorBase<VectorType>::value_type
+BlockVectorBase<VectorType>::
+add_and_dot (const typename BlockVectorBase<VectorType>::value_type a,
+             const BlockVectorBase<VectorType> &V,
+             const BlockVectorBase<VectorType> &W)
+{
+  AssertDimension (n_blocks(), V.n_blocks());
+  AssertDimension (n_blocks(), W.n_blocks());
+
+  value_type sum = 0.;
+  for (size_type i=0; i<n_blocks(); ++i)
+    sum += components[i].add_and_dot(a, V.components[i], W.components[i]);
+
+  return sum;
+}
+
+
+
 template <class VectorType>
 BlockVectorBase<VectorType> &
 BlockVectorBase<VectorType>::operator += (const BlockVectorBase<VectorType> &v)
index 2c38416dfb0640b83f54fa7c2b1efbad773e8b92..2ca70988143b73c8ce06aa6d11138b8baf9ffb94 100644 (file)
@@ -377,6 +377,26 @@ namespace parallel
        */
       real_type linfty_norm () const;
 
+      /**
+       * 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
+       *
+       * The reason this function exists is that this operation involves less
+       * memory transfer than calling the two functions separately. This
+       * method only needs to load three vectors, @p this, @p V, @p W, whereas
+       * calling separate methods means to load the calling vector @p this
+       * twice. Since most vector operations are memory transfer limited, this
+       * reduces the time by 25\% (or 50\% if @p W equals @p this).
+       */
+      Number add_and_dot (const Number               a,
+                          const BlockVector<Number> &V,
+                          const BlockVector<Number> &W);
+
       /**
        * Scale each element of the vector by the given factor.
        *
@@ -928,7 +948,28 @@ namespace parallel
 
     template <typename Number>
     inline
-    void BlockVector<Number>::swap (BlockVector<Number> &v)
+    Number
+    BlockVector<Number>::add_and_dot (const Number               a,
+                                      const BlockVector<Number> &V,
+                                      const BlockVector<Number> &W)
+    {
+      Number local_result = Number();
+      for (unsigned int i=0; i<this->n_blocks(); ++i)
+        local_result += this->block(i).add_and_dot_local(a, V.block(i), W.block(i));
+
+      if (this->block(0).partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    this->block(0).partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
+    template <typename Number>
+    inline
+    void
+    BlockVector<Number>::swap (BlockVector<Number> &v)
     {
       Assert (this->n_blocks() == v.n_blocks(),
               ExcDimensionMismatch(this->n_blocks(), v.n_blocks()));
index 416d615f1d5464d6e4937f732490aef684b3c913..06c64f01ece4dc6b2594812729911c163f26a14b 100644 (file)
@@ -561,6 +561,26 @@ namespace parallel
        */
       real_type linfty_norm () const;
 
+      /**
+       * 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
+       *
+       * The reason this function exists is that this operation involves less
+       * memory transfer than calling the two functions separately. This
+       * method only needs to load three vectors, @p this, @p V, @p W, whereas
+       * calling separate methods means to load the calling vector @p this
+       * twice. Since most vector operations are memory transfer limited, this
+       * reduces the time by 25\% (or 50\% if @p W equals @p this).
+       */
+      Number add_and_dot (const Number          a,
+                          const Vector<Number> &V,
+                          const Vector<Number> &W);
+
       /**
        * Returns the global size of the vector, equal to the sum of the number
        * of locally owned indices among all the processors.
@@ -1009,6 +1029,14 @@ namespace parallel
        */
       real_type linfty_norm_local () const;
 
+      /**
+       * Local part of the addition followed by an inner product of two
+       * vectors.
+       */
+      Number add_and_dot_local (const Number          a,
+                                const Vector<Number> &V,
+                                const Vector<Number> &W);
+
       /**
        * Shared pointer to store the parallel partitioning information. This
        * information can be shared between several vectors that have the same
@@ -1203,9 +1231,7 @@ namespace parallel
     inline
     Vector<Number>::~Vector ()
     {
-      if (val != 0)
-        delete[] val;
-      val = 0;
+      resize_val(0);
 
       if (import_data != 0)
         delete[] import_data;
@@ -1630,6 +1656,39 @@ namespace parallel
 
 
 
+    template <typename Number>
+    inline
+    Number
+    Vector<Number>::add_and_dot_local(const Number          a,
+                                      const Vector<Number> &V,
+                                      const Vector<Number> &W)
+    {
+      // on some processors, the size might be zero, which is not allowed by
+      // the dealii::Vector class. Therefore, insert a check here
+      return (partitioner->local_size()>0 ?
+              vector_view.add_and_dot(a, V.vector_view, W.vector_view)
+              : Number());
+    }
+
+
+
+    template <typename Number>
+    inline
+    Number
+    Vector<Number>::add_and_dot (const Number          a,
+                                 const Vector<Number> &V,
+                                 const Vector<Number> &W)
+    {
+      Number local_result = add_and_dot_local(a, V, W);
+      if (partitioner->n_mpi_processes() > 1)
+        return Utilities::MPI::sum (local_result,
+                                    partitioner->get_communicator());
+      else
+        return local_result;
+    }
+
+
+
     template <typename Number>
     inline
     typename Vector<Number>::size_type
index 6b2c55422247c4a80c98bee6d8447f1a15a0db54..bd5b84af25b84e143df9d4fc0f9749ce270092ec 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 2011 - 2013 by the deal.II authors
+// Copyright (C) 2011 - 2014 by the deal.II authors
 //
 // This file is part of the deal.II library.
 //
@@ -24,6 +24,8 @@
 #include <deal.II/lac/petsc_parallel_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 
+#include <mm_malloc.h>
+
 DEAL_II_NAMESPACE_OPEN
 
 
@@ -57,14 +59,14 @@ namespace parallel
           Assert (((allocated_size > 0 && val != 0) ||
                    val == 0), ExcInternalError());
           if (val != 0)
-            delete [] val;
-          val = new Number[new_alloc_size];
+            _mm_free(val);
+          val = static_cast<Number *>(_mm_malloc (sizeof(Number)*new_alloc_size, 64));
           allocated_size = new_alloc_size;
         }
       else if (new_alloc_size == 0)
         {
           if (val != 0)
-            delete [] val;
+            _mm_free(val);
           val = 0;
           allocated_size = 0;
         }
index efb6b08201a61a6eae187dbbfee48f17a00420eb..8a67bf503e79f13dcb0055bd9572ed9ee4bd184c 100644 (file)
@@ -572,6 +572,25 @@ namespace PETScWrappers
      */
     real_type linfty_norm () const;
 
+    /**
+     * 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
+     *
+     * The reason this function exists is for compatibility with deal.II's own
+     * vector classes which can implement this functionality with less memory
+     * transfer. However, for PETSc vectors such a combined operation is not
+     * natively supported and thus the cost is completely equivalent as
+     * calling the two methods separately.
+     */
+    PetscScalar add_and_dot (const Number      a,
+                             const VectorBase &V,
+                             const VectorBase &W);
+
     /**
      * Normalize vector by dividing by the $l_2$-norm of the
      * vector. Return the vector norm before normalization.
index 5041c67c19e3a3f03dba2acd9bd391aa672c0aed..f86d60e47e480b2a6f432c50217f16e424a6ab94 100644 (file)
@@ -351,11 +351,11 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX &A,
       if (std::fabs(alpha) > 1.e10)
         return true;
 
-      r.add(-alpha, v);
+      res = std::sqrt(r.add_and_dot(-alpha, v, r));
 
       // check for early success, see the lac/bicgstab_early testcase as to
       // why this is necessary
-      if (this->control().check(step, r.l2_norm()) == SolverControl::success)
+      if (this->control().check(step, res) == SolverControl::success)
         {
           Vx->add(alpha, y);
           print_vectors(step, *Vx, r, y);
@@ -367,12 +367,14 @@ SolverBicgstab<VECTOR>::iterate(const MATRIX &A,
       rhobar = t*r;
       omega = rhobar/(t*t);
       Vx->add(alpha, y, omega, z);
-      r.add(-omega, t);
 
       if (additional_data.exact_residual)
-        res = criterion(A, *Vx, *Vb);
+        {
+          r.add(-omega, t);
+          res = criterion(A, *Vx, *Vb);
+        }
       else
-        res = r.l2_norm();
+        res = std::sqrt(r.add_and_dot(-omega, t, r));
 
       state = this->control().check(step, res);
       print_vectors(step, *Vx, r, y);
index 9f57f281278d01793d856f7726f28fd3a4000289..9facc394eacffa24d5adbe7930e8ee6c24c590d9 100644 (file)
@@ -391,9 +391,8 @@ SolverCG<VECTOR>::solve (const MATRIX         &A,
           Assert(alpha != 0., ExcDivideByZero());
           alpha = gh/alpha;
 
-          g.add(alpha,h);
           x.add(alpha,d);
-          res = g.l2_norm();
+          res = std::sqrt(g.add_and_dot(alpha, h, g));
 
           print_vectors(it, x, g, d);
 
index 9e4d9333258b5cafe563c54c15051918c6ff065a..a5c0622ed653281f7421cc8d792e34423af39d6f 100644 (file)
@@ -528,6 +528,7 @@ SolverGMRES<VECTOR>::modified_gram_schmidt (const internal::SolverGMRES::TmpVect
                                             Vector<double>     &h,
                                             bool               &re_orthogonalize)
 {
+  Assert(dim > 0, ExcInternalError());
   const unsigned int inner_iteration = dim - 1;
 
   // need initial norm for detection of re-orthogonalization, see below
@@ -536,11 +537,10 @@ SolverGMRES<VECTOR>::modified_gram_schmidt (const internal::SolverGMRES::TmpVect
     norm_vv_start = vv.l2_norm();
 
   // Orthogonalization
-  for (unsigned int i=0 ; i<dim ; ++i)
-    {
-      h(i) = vv * orthogonal_vectors[i];
-      vv.add(-h(i), orthogonal_vectors[i]);
-    };
+  h(0) = vv * orthogonal_vectors[0];
+  for (unsigned int i=1 ; i<dim ; ++i)
+    h(i) = vv.add_and_dot(-h(i-1), orthogonal_vectors[i-1], orthogonal_vectors[i]);
+  double norm_vv = std::sqrt(vv.add_and_dot(-h(dim-1), orthogonal_vectors[dim-1], vv));
 
   // Re-orthogonalization if loss of orthogonality detected. For the test, use
   // a strategy discussed in C. T. Kelley, Iterative Methods for Linear and
@@ -551,7 +551,6 @@ SolverGMRES<VECTOR>::modified_gram_schmidt (const internal::SolverGMRES::TmpVect
   // previous vectors, which indicates loss of precision.
   if (re_orthogonalize == false && inner_iteration % 5 == 4)
     {
-      const double norm_vv = vv.l2_norm();
       if (norm_vv > 10. * norm_vv_start *
           std::sqrt(std::numeric_limits<typename VECTOR::value_type>::epsilon()))
         return norm_vv;
@@ -565,14 +564,18 @@ SolverGMRES<VECTOR>::modified_gram_schmidt (const internal::SolverGMRES::TmpVect
     }
 
   if (re_orthogonalize == true)
-    for (unsigned int i=0 ; i<dim ; ++i)
-      {
-        double htmp = vv * orthogonal_vectors[i];
-        h(i) += htmp;
-        vv.add(-htmp, orthogonal_vectors[i]);
-      }
-
-  return vv.l2_norm();
+    {
+      double htmp = vv * orthogonal_vectors[0];
+      h(0) += htmp;
+      for (unsigned int i=1 ; i<dim ; ++i)
+        {
+          htmp = vv.add_and_dot(-htmp, orthogonal_vectors[i-1], orthogonal_vectors[i]);
+          h(i) += htmp;
+        }
+      norm_vv = std::sqrt(vv.add_and_dot(-htmp, orthogonal_vectors[dim-1], vv));
+    }
+
+  return norm_vv;
 }
 
 
@@ -980,12 +983,10 @@ SolverFGMRES<VECTOR>::solve (
           A.vmult(*aux, z[j]);
 
           // Gram-Schmidt
-          for (unsigned int i=0; i<=j; ++i)
-            {
-              H(i,j) = *aux * v[i];
-              aux->add(-H(i,j), v[i]);
-            }
-          H(j+1,j) = a = aux->l2_norm();
+          H(0,j) = *aux * v[0];
+          for (unsigned int i=1; i<=j; ++i)
+            H(i,j) = aux->add_and_dot(-H(i-1,j), v[i-1], v[i]);
+          H(j+1,j) = a = std::sqrt(aux->add_and_dot(-H(j,j), v[j], *aux));
 
           // Compute projected solution
 
index 449c4a5cf05cc573d3c586b229049fe369512601..5fbdf728ec305168c9bbdebaa621fe8ea2953e1e 100644 (file)
@@ -468,6 +468,25 @@ namespace TrilinosWrappers
      */
     real_type linfty_norm () const;
 
+    /**
+     * 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
+     *
+     * The reason this function exists is for compatibility with deal.II's own
+     * vector classes which can implement this functionality with less memory
+     * transfer. However, for Trilinos vectors such a combined operation is
+     * not natively supported and thus the cost is completely equivalent as
+     * calling the two methods separately.
+     */
+    TrilinosScalar add_and_dot (const TrilinosScalar a,
+                                const VectorBase    &V,
+                                const VectorBase    &W);
+
     /**
      * Return whether the vector contains only elements with value
      * zero. This is a collective operation. This function is expensive, because
@@ -1565,6 +1584,18 @@ namespace TrilinosWrappers
 
 
 
+  inline
+  TrilinosScalar
+  VectorBase::add_and_dot (const TrilinosScalar a,
+                           const VectorBase &V,
+                           const VectorBase &W)
+  {
+    this->add(a, V);
+    return *this * W;
+  }
+
+
+
   // inline also scalar products, vector
   // additions etc. since they are all
   // representable by a single Trilinos
index 47765ee3dc00907dc5beae767238894df087dab4..4c037a8e3b59b2008e2fae54115d21f86d3c6918 100644 (file)
@@ -489,6 +489,31 @@ public:
    * Maximum absolute value of the elements.
    */
   real_type linfty_norm () const;
+
+  /**
+   * 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
+   *
+   * The reason this function exists is that this operation involves less
+   * memory transfer than calling the two functions separately. This method
+   * only needs to load three vectors, @p this, @p V, @p W, whereas calling
+   * separate methods means to load the calling vector @p this twice. Since
+   * most vector operations are memory transfer limited, this reduces the time
+   * by 25\% (or 50\% if @p W equals @p this).
+   *
+   * @dealiiOperationIsMultithreaded The algorithm uses pairwise
+   * summation with the same order of summation in every run, which gives
+   * fully repeatable results from one run to another.
+   */
+  Number add_and_dot (const Number          a,
+                      const Vector<Number> &V,
+                      const Vector<Number> &W);
+
   //@}
 
 
index d6073f503080348724e634a794aedfbfa1a739d5..a331a19a8385aa526b28621b9babd2ef2513df7e 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/base/parallel.h>
 #include <deal.II/base/thread_management.h>
 #include <deal.II/base/multithread_info.h>
+#include <deal.II/base/vectorization.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/block_vector.h>
 
@@ -823,57 +824,224 @@ namespace internal
   namespace Vector
   {
     // All sums over all the vector entries (l2-norm, inner product, etc.) are
-    // performed with the same code, using a templated operation defined here
+    // performed with the same code, using a templated operation defined
+    // here. There are always two versions defined, a standard one that covers
+    // most cases and a vectorized one which is only for equal types and float
+    // and double.
     template <typename Number, typename Number2>
-    struct InnerProd
+    struct Dot
     {
       Number
-      operator() (const Number *&X, const Number2 *&Y, const Number &) const
+      operator() (const Number *&X, const Number2 *&Y, const Number &, Number *&) const
       {
         return *X++ * Number(numbers::NumberTraits<Number2>::conjugate(*Y++));
       }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&X, const Number *&Y, const Number &, Number *&) const
+      {
+        VectorizedArray<Number> x, y;
+        x.load(X);
+        y.load(Y);
+        X += VectorizedArray<Number>::n_array_elements;
+        Y += VectorizedArray<Number>::n_array_elements;
+        return x * y;
+      }
     };
 
     template <typename Number, typename RealType>
     struct Norm2
     {
       RealType
-      operator() (const Number  *&X, const Number  *&, const RealType &) const
+      operator() (const Number  *&X, const Number  *&, const RealType &, Number *&) const
       {
         return numbers::NumberTraits<Number>::abs_square(*X++);
       }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const
+      {
+        VectorizedArray<Number> x;
+        x.load(X);
+        X += VectorizedArray<Number>::n_array_elements;
+        return x * x;
+      }
     };
 
     template <typename Number, typename RealType>
     struct Norm1
     {
       RealType
-      operator() (const Number  *&X, const Number  *&, const RealType &) const
+      operator() (const Number  *&X, const Number  *&, const RealType &, Number *&) const
       {
         return numbers::NumberTraits<Number>::abs(*X++);
       }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const
+      {
+        VectorizedArray<Number> x;
+        x.load(X);
+        X += VectorizedArray<Number>::n_array_elements;
+        return std::abs(x);
+      }
     };
 
     template <typename Number, typename RealType>
     struct NormP
     {
       RealType
-      operator() (const Number  *&X, const Number  *&, const RealType &p) const
+      operator() (const Number  *&X, const Number  *&, const RealType &p, Number *&) const
       {
         return std::pow(numbers::NumberTraits<Number>::abs(*X++), p);
       }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&X, const Number *&, const Number &p, Number *&) const
+      {
+        VectorizedArray<Number> x;
+        x.load(X);
+        X += VectorizedArray<Number>::n_array_elements;
+        return std::pow(std::abs(x),p);
+      }
     };
 
     template <typename Number>
     struct MeanValue
     {
       Number
-      operator() (const Number  *&X, const Number  *&, const Number &) const
+      operator() (const Number  *&X, const Number  *&, const Number &, Number *&) const
       {
         return *X++;
       }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const
+      {
+        VectorizedArray<Number> x;
+        x.load(X);
+        X += VectorizedArray<Number>::n_array_elements;
+        return x;
+      }
     };
 
+    template <typename Number>
+    struct AddAndDot
+    {
+      Number
+      operator() (const Number *&V, const Number *&W, const Number &a,
+                  Number *&X) const
+      {
+        *X += a **V++;
+        return *X++ * Number(numbers::NumberTraits<Number>::conjugate(*W++));
+      }
+
+      VectorizedArray<Number>
+      do_vectorized(const Number *&V, const Number *&W, const Number &a,
+                    Number *&X) const
+      {
+        VectorizedArray<Number> x, w, v;
+        x.load(X);
+        v.load(V);
+        x += a * v;
+        x.store(X);
+        // may only load from W after storing in X because the pointers might
+        // point to the same memory
+        w.load(W);
+        X += VectorizedArray<Number>::n_array_elements;
+        V += VectorizedArray<Number>::n_array_elements;
+        W += VectorizedArray<Number>::n_array_elements;
+        return x * w;
+      }
+    };
+
+
+
+    // this is the inner working routine for the accumulation loops
+    // below. This is the standard case where the loop bounds are known. We
+    // pulled this function out of the regular accumulate routine because we
+    // might do this thing vectorized (see specialized function below)
+    template <typename Operation, typename Number, typename Number2,
+              typename ResultType, typename size_type>
+    void
+    accumulate_regular(const Operation &op,
+                       const Number   *&X,
+                       const Number2  *&Y,
+                       const ResultType power,
+                       const size_type  n_chunks,
+                       Number         *&Z,
+                       ResultType (&outer_results)[128],
+                       internal::bool2type<false>,
+                       const unsigned int start_chunk=0)
+    {
+      AssertIndexRange(start_chunk, n_chunks+1);
+      for (size_type i=start_chunk; i<n_chunks; ++i)
+        {
+          ResultType r0 = op(X, Y, power, Z);
+          ResultType r1 = op(X, Y, power, Z);
+          ResultType r2 = op(X, Y, power, Z);
+          ResultType r3 = op(X, Y, power, Z);
+          for (size_type j=1; j<8; ++j)
+            {
+              r0 += op(X, Y, power, Z);
+              r1 += op(X, Y, power, Z);
+              r2 += op(X, Y, power, Z);
+              r3 += op(X, Y, power, Z);
+            }
+          r0 += r1;
+          r2 += r3;
+          outer_results[i] = r0 + r2;
+        }
+    }
+
+
+
+    // this is the inner working routine for the accumulation loops
+    // below. This is the specialized case where the loop bounds are known and
+    // where we can vectorize. In that case, we request the 'do_vectorized'
+    // routine of the operation instead of the regular one which does several
+    // operations at once.
+    template <typename Operation, typename Number, typename size_type>
+    void
+    accumulate_regular(const Operation &op,
+                       const Number   *&X,
+                       const Number   *&Y,
+                       const Number     power,
+                       const size_type  n_chunks,
+                       Number         *&Z,
+                       Number (&outer_results)[128],
+                       internal::bool2type<true>)
+    {
+      for (size_type i=0; i<n_chunks/VectorizedArray<Number>::n_array_elements; ++i)
+        {
+          VectorizedArray<Number> r0 = op.do_vectorized(X, Y, power, Z);
+          VectorizedArray<Number> r1 = op.do_vectorized(X, Y, power, Z);
+          VectorizedArray<Number> r2 = op.do_vectorized(X, Y, power, Z);
+          VectorizedArray<Number> r3 = op.do_vectorized(X, Y, power, Z);
+          for (size_type j=1; j<8; ++j)
+            {
+              r0 += op.do_vectorized(X, Y, power, Z);
+              r1 += op.do_vectorized(X, Y, power, Z);
+              r2 += op.do_vectorized(X, Y, power, Z);
+              r3 += op.do_vectorized(X, Y, power, Z);
+            }
+          r0 += r1;
+          r2 += r3;
+          r0 += r2;
+          r0.store(&outer_results[i*VectorizedArray<Number>::n_array_elements]);
+        }
+
+      // If we are treating a case where the vector length is not divisible by
+      // the vectorization length, need the other routine for the cleanup work
+      if (n_chunks % VectorizedArray<Number>::n_array_elements != 0)
+        accumulate_regular(op, X, Y, power, n_chunks, Z, outer_results,
+                           internal::bool2type<false>(),
+                           (n_chunks/VectorizedArray<Number>::n_array_elements)
+                           * VectorizedArray<Number>::n_array_elements);
+    }
+
+
+
     // this is the main working loop for all vector sums using the templated
     // operation above. it accumulates the sums using a block-wise summation
     // algorithm with post-update. this blocked algorithm has been proposed in
@@ -903,6 +1071,7 @@ namespace internal
                      const Number2     *Y,
                      const ResultType   power,
                      const size_type    vec_size,
+                     Number            *Z,
                      ResultType        &result,
                      const int          depth = -1)
     {
@@ -917,29 +1086,17 @@ namespace internal
           const size_type remainder = vec_size % 32;
           Assert (remainder == 0 || n_chunks < 128, ExcInternalError());
 
-          for (size_type i=0; i<n_chunks; ++i)
-            {
-              ResultType r0 = op(X, Y, power);
-              for (size_type j=1; j<8; ++j)
-                r0 += op(X, Y, power);
-              ResultType r1 = op(X, Y, power);
-              for (size_type j=1; j<8; ++j)
-                r1 += op(X, Y, power);
-              r0 += r1;
-              r1 = op(X, Y, power);
-              for (size_type j=1; j<8; ++j)
-                r1 += op(X, Y, power);
-              ResultType r2 = op(X, Y, power);
-              for (size_type j=1; j<8; ++j)
-                r2 += op(X, Y, power);
-              r1 += r2;
-              r0 += r1;
-              outer_results[i] = r0;
-            }
-
-          // now work on the remainder, i.e., the last
-          // up to 32 values. Use switch statement with
-          // fall-through to work on these values.
+          // Select between the regular version and vectorized version based
+          // on the number types we are given. To choose the vectorized
+          // version often enough, we need to have all tasks but the last one
+          // to be divisible by the vectorization length
+          accumulate_regular(op, X, Y, power, n_chunks, Z, outer_results,
+                             internal::bool2type<(types_are_equal<Number,Number2>::value &&
+                                                  (types_are_equal<Number,double>::value ||
+                                                   types_are_equal<Number,float>::value))>());
+
+          // now work on the remainder, i.e., the last up to 32 values. Use
+          // switch statement with fall-through to work on these values.
           if (remainder > 0)
             {
               const size_type inner_chunks = remainder / 8;
@@ -950,24 +1107,24 @@ namespace internal
               switch (inner_chunks)
                 {
                 case 3:
-                  r2 = op(X, Y, power);
+                  r2 = op(X, Y, power, Z);
                   for (size_type j=1; j<8; ++j)
-                    r2 += op(X, Y, power);
+                    r2 += op(X, Y, power, Z);
                 // no break
                 case 2:
-                  r1 = op(X, Y, power);
+                  r1 = op(X, Y, power, Z);
                   for (size_type j=1; j<8; ++j)
-                    r1 += op(X, Y, power);
+                    r1 += op(X, Y, power, Z);
                   r1 += r2;
                 // no break
                 case 1:
-                  r2 = op(X, Y, power);
+                  r2 = op(X, Y, power, Z);
                   for (size_type j=1; j<8; ++j)
-                    r2 += op(X, Y, power);
+                    r2 += op(X, Y, power, Z);
                 // no break
                 default:
                   for (size_type j=0; j<remainder_inner; ++j)
-                    r0 += op(X, Y, power);
+                    r0 += op(X, Y, power, Z);
                   r0 += r2;
                   r0 += r1;
                   outer_results[n_chunks] = r0;
@@ -1010,19 +1167,20 @@ namespace internal
           Threads::TaskGroup<> task_group;
           task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
                                           ResultType,size_type>,
-                                          op, X, Y, power, new_size, r0, next_depth);
+                                          op, X, Y, power, new_size, Z, r0, next_depth);
           task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
                                           ResultType,size_type>,
                                           op, X+new_size, Y+new_size, power,
-                                          new_size, r1, next_depth);
+                                          new_size, Z+new_size, r1, next_depth);
           task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
                                           ResultType,size_type>,
                                           op, X+2*new_size, Y+2*new_size, power,
-                                          new_size, r2, next_depth);
+                                          new_size, Z+2*new_size, r2, next_depth);
           task_group += Threads::new_task(&accumulate<Operation,Number,Number2,
                                           ResultType,size_type>,
                                           op, X+3*new_size, Y+3*new_size, power,
-                                          vec_size-3*new_size, r3, next_depth);
+                                          vec_size-3*new_size, Z+3*new_size, r3,
+                                          next_depth);
           task_group.join_all();
           r0 += r1;
           r2 += r3;
@@ -1031,15 +1189,15 @@ namespace internal
 #endif
       else
         {
-          // split vector into four pieces and work on
-          // the pieces recursively. Make pieces (except last)
-          // divisible by 1024.
+          // split vector into four pieces and work on the pieces
+          // recursively. Make pieces (except last) divisible by 1024.
           const size_type new_size = (vec_size / 4096) * 1024;
           ResultType r0, r1, r2, r3;
-          accumulate (op, X, Y, power, new_size, r0);
-          accumulate (op, X+new_size, Y+new_size, power, new_size, r1);
-          accumulate (op, X+2*new_size, Y+2*new_size, power, new_size, r2);
-          accumulate (op, X+3*new_size, Y+3*new_size, power, vec_size-3*new_size, r3);
+          accumulate (op, X, Y, power, new_size, Z, r0);
+          accumulate (op, X+new_size, Y+new_size, power, new_size, Z+new_size, r1);
+          accumulate (op, X+2*new_size, Y+2*new_size, power, new_size, Z+2*new_size, r2);
+          accumulate (op, X+3*new_size, Y+3*new_size, power, vec_size-3*new_size,
+                      Z+3*new_size, r3);
           r0 += r1;
           r2 += r3;
           result = r0 + r2;
@@ -1063,14 +1221,15 @@ Number Vector<Number>::operator * (const Vector<Number2> &v) const
           ExcDimensionMismatch(vec_size, v.size()));
 
   Number sum;
-  internal::Vector::accumulate (internal::Vector::InnerProd<Number,Number2>(),
-                                val, v.val, Number(), vec_size, sum);
+  internal::Vector::accumulate (internal::Vector::Dot<Number,Number2>(),
+                                val, v.val, Number(), vec_size, val, sum);
   Assert(numbers::is_finite(sum), ExcNumberNotFinite());
 
   return sum;
 }
 
 
+
 template <typename Number>
 typename Vector<Number>::real_type
 Vector<Number>::norm_sqr () const
@@ -1079,7 +1238,7 @@ Vector<Number>::norm_sqr () const
 
   real_type sum;
   internal::Vector::accumulate (internal::Vector::Norm2<Number,real_type>(),
-                                val, val, real_type(), vec_size, sum);
+                                val, val, real_type(), vec_size, val, sum);
 
   Assert(numbers::is_finite(sum), ExcNumberNotFinite());
 
@@ -1087,6 +1246,7 @@ Vector<Number>::norm_sqr () const
 }
 
 
+
 template <typename Number>
 Number Vector<Number>::mean_value () const
 {
@@ -1094,7 +1254,7 @@ Number Vector<Number>::mean_value () const
 
   Number sum;
   internal::Vector::accumulate (internal::Vector::MeanValue<Number>(),
-                                val, val, Number(), vec_size, sum);
+                                val, val, Number(), vec_size, val, sum);
 
   return sum / real_type(size());
 }
@@ -1109,7 +1269,7 @@ Vector<Number>::l1_norm () const
 
   real_type sum;
   internal::Vector::accumulate (internal::Vector::Norm1<Number,real_type>(),
-                                val, val, real_type(), vec_size, sum);
+                                val, val, real_type(), vec_size, val, sum);
 
   return sum;
 }
@@ -1129,7 +1289,7 @@ Vector<Number>::l2_norm () const
 
   real_type norm_square;
   internal::Vector::accumulate (internal::Vector::Norm2<Number,real_type>(),
-                                val, val, real_type(), vec_size, norm_square);
+                                val, val, real_type(), vec_size, val, norm_square);
   if (numbers::is_finite(norm_square) &&
       norm_square >= std::numeric_limits<real_type>::min())
     return std::sqrt(norm_square);
@@ -1172,7 +1332,7 @@ Vector<Number>::lp_norm (const real_type p) const
 
   real_type sum;
   internal::Vector::accumulate (internal::Vector::NormP<Number,real_type>(),
-                                val, val, p, vec_size, sum);
+                                val, val, p, vec_size, val, sum);
 
   if (numbers::is_finite(sum) && sum >= std::numeric_limits<real_type>::min())
     return std::pow(sum, static_cast<real_type>(1./p));
@@ -1216,6 +1376,27 @@ Vector<Number>::linfty_norm () const
 }
 
 
+
+template <typename Number>
+Number
+Vector<Number>::add_and_dot (const Number          a,
+                             const Vector<Number> &V,
+                             const Vector<Number> &W)
+{
+  Assert (vec_size!=0, ExcEmptyObject());
+  AssertDimension (vec_size, V.size());
+  AssertDimension (vec_size, W.size());
+
+  Number sum;
+  internal::Vector::accumulate (internal::Vector::AddAndDot<Number>(),
+                                V.val, W.val, a, vec_size, val, sum);
+  Assert(numbers::is_finite(sum), ExcNumberNotFinite());
+
+  return sum;
+}
+
+
+
 template <typename Number>
 Vector<Number> &Vector<Number>::operator += (const Vector<Number> &v)
 {
@@ -1226,6 +1407,7 @@ Vector<Number> &Vector<Number>::operator += (const Vector<Number> &v)
 }
 
 
+
 template <typename Number>
 Vector<Number> &Vector<Number>::operator -= (const Vector<Number> &v)
 {
@@ -1241,6 +1423,7 @@ Vector<Number> &Vector<Number>::operator -= (const Vector<Number> &v)
 }
 
 
+
 template <typename Number>
 void Vector<Number>::add (const Number v)
 {
@@ -1253,6 +1436,7 @@ void Vector<Number>::add (const Number v)
 }
 
 
+
 template <typename Number>
 void Vector<Number>::add (const Vector<Number> &v)
 {
@@ -1266,6 +1450,7 @@ void Vector<Number>::add (const Vector<Number> &v)
 }
 
 
+
 template <typename Number>
 void Vector<Number>::add (const Number a, const Vector<Number> &v,
                           const Number b, const Vector<Number> &w)
index 041648c52e74ed490f1a669cd78532c2b851c3f9..26324e95f69ebc56c6ec0b4f61460f45fa2dbbff 100644 (file)
@@ -376,6 +376,17 @@ namespace PETScWrappers
 
 
 
+  PetscScalar
+  VectorBase::add_and_dot (const PetscScalar a,
+                           const VectorBase &V,
+                           const VectorBase &W)
+  {
+    this->add(a, V);
+    return *this * W;
+  }
+
+
+
   void
   VectorBase::compress (::dealii::VectorOperation::values)
   {
@@ -825,7 +836,7 @@ namespace PETScWrappers
 
   void
   VectorBase::add (const PetscScalar a,
-                   const VectorBase     &v)
+                   const VectorBase &v)
   {
     Assert (!has_ghost_elements(), ExcGhostsPresent());
     Assert (numbers::is_finite(a), ExcNumberNotFinite());
index f209d1462e0ce1fc2976bd643677549c3dca046e..bb310136036d48f55c5bdb8669fbc5cf19ab7cfc 100644 (file)
@@ -40,7 +40,7 @@ void test ()
   for (unsigned int i=0; i<n; ++i)
     matrix.diag_element(i) = (i+1);
 
-  SolverControl control(1000, 1e2*std::numeric_limits<number>::epsilon());
+  SolverControl control(1000, 1e3*std::numeric_limits<number>::epsilon());
   typename SolverGMRES<Vector<number> >::AdditionalData data;
   data.max_n_tmp_vectors = 202;
 
@@ -59,7 +59,7 @@ int main()
   deallog.push("double");
   test<double>();
   deallog.pop();
-  deallog.threshold_double(1.e-4);
+  deallog.threshold_double(2.e-4);
   deallog.push("float");
   test<float>();
   deallog.pop();
diff --git a/tests/lac/gmres_reorthogonalize_02.compiler=Intel=no.release.output b/tests/lac/gmres_reorthogonalize_02.compiler=Intel=no.release.output
deleted file mode 100644 (file)
index 561b5e4..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL:double:GMRES::Starting value 14.14
-DEAL:double:GMRES::Convergence step 109 value 0
-DEAL:float:GMRES::Starting value 14.14
-DEAL:float:GMRES::Convergence step 68 value 0
diff --git a/tests/lac/gmres_reorthogonalize_02.compiler=Intel=yes.release.output b/tests/lac/gmres_reorthogonalize_02.compiler=Intel=yes.release.output
deleted file mode 100644 (file)
index e187d5b..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL:double:GMRES::Starting value 14.14
-DEAL:double:GMRES::Convergence step 109 value 0
-DEAL:float:GMRES::Starting value 14.14
-DEAL:float:GMRES::Convergence step 69 value 0
diff --git a/tests/lac/gmres_reorthogonalize_02.debug.output b/tests/lac/gmres_reorthogonalize_02.debug.output
deleted file mode 100644 (file)
index 561b5e4..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL:double:GMRES::Starting value 14.14
-DEAL:double:GMRES::Convergence step 109 value 0
-DEAL:float:GMRES::Starting value 14.14
-DEAL:float:GMRES::Convergence step 68 value 0
diff --git a/tests/lac/gmres_reorthogonalize_02.output b/tests/lac/gmres_reorthogonalize_02.output
new file mode 100644 (file)
index 0000000..9b62b4f
--- /dev/null
@@ -0,0 +1,5 @@
+
+DEAL:double:GMRES::Starting value 14.14
+DEAL:double:GMRES::Convergence step 105 value 0
+DEAL:float:GMRES::Starting value 14.14
+DEAL:float:GMRES::Convergence step 59 value 0
index b021e983feaf9530eac7e3db1bd7e729f10efdfe..344671676a63cf98c7e947403c845d3dd814d36d 100644 (file)
@@ -87,7 +87,7 @@ int main()
   deallog.push("double");
   test<double>();
   deallog.pop();
-  deallog.threshold_double(1.e-4);
+  deallog.threshold_double(2.e-4);
   deallog.push("float");
   test<float>();
   deallog.pop();
diff --git a/tests/lac/gmres_reorthogonalize_04.compiler=Intel=no.release.output b/tests/lac/gmres_reorthogonalize_04.compiler=Intel=no.release.output
deleted file mode 100644 (file)
index 06bb180..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL:double:GMRES::Starting value 14.1
-DEAL:double:GMRES::Convergence step 105 value 0
-DEAL:float:GMRES::Starting value 14.1
-DEAL:float:GMRES::Convergence step 59 value 0.000114
diff --git a/tests/lac/gmres_reorthogonalize_04.debug.output b/tests/lac/gmres_reorthogonalize_04.debug.output
deleted file mode 100644 (file)
index 06bb180..0000000
+++ /dev/null
@@ -1,5 +0,0 @@
-
-DEAL:double:GMRES::Starting value 14.1
-DEAL:double:GMRES::Convergence step 105 value 0
-DEAL:float:GMRES::Starting value 14.1
-DEAL:float:GMRES::Convergence step 59 value 0.000114
similarity index 70%
rename from tests/lac/gmres_reorthogonalize_04.compiler=Intel=yes.release.output
rename to tests/lac/gmres_reorthogonalize_04.output
index 84f47a47656eab87da024d540693b12a06319b69..5d6219e85115156c9484fe8463cf9ae81549bdd4 100644 (file)
@@ -2,4 +2,4 @@
 DEAL:double:GMRES::Starting value 14.1
 DEAL:double:GMRES::Convergence step 105 value 0
 DEAL:float:GMRES::Starting value 14.1
-DEAL:float:GMRES::Convergence step 59 value 0.000110
+DEAL:float:GMRES::Convergence step 59 value 0
diff --git a/tests/lac/vector_accumulation_02.cc b/tests/lac/vector_accumulation_02.cc
new file mode 100644 (file)
index 0000000..9dfa755
--- /dev/null
@@ -0,0 +1,82 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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 the inner product is exactly the same down to roundoff for
+// various vectors. Due to vectorization and its requirements for alignment,
+// we also test that the result does not depend on the alignment of any of the
+// vectors by choosing among several start locations of a vector relative to a
+// data array allocated as usual.
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/vector_view.h>
+#include <cmath>
+#include <fstream>
+#include <iomanip>
+
+
+
+
+template <typename number>
+void check_norms ()
+{
+  for (unsigned int test=0; test<5; ++test)
+    {
+      const unsigned int size = Testing::rand() % 20000;
+      Vector<number> larger1 (size+8), larger2(size+8), in1(size), in2(size);
+      for (unsigned int i=0; i<size; ++i)
+        {
+          in1(i) = static_cast<number>(Testing::rand())/static_cast<number>(RAND_MAX);
+          in2(i) = static_cast<number>(Testing::rand())/static_cast<number>(RAND_MAX);
+        }
+
+      const number inner_product = in1 * in2;
+      deallog << "Vector difference: ";
+      for (unsigned int shift1=0; shift1<8; ++shift1)
+        for (unsigned int shift2=0; shift2<8; ++shift2)
+          {
+            VectorView<number> v1 (size, larger1.begin()+shift1);
+            VectorView<number> v2 (size, larger2.begin()+shift2);
+            for (unsigned int i=0; i<size; ++i)
+              {
+                v1(i) = in1(i);
+                v2(i) = in2(i);
+              }
+
+            const number result = v1 * v2;
+            deallog << static_cast<double>(std::abs(result-inner_product)) << " ";
+          }
+      deallog << std::endl;
+    }
+}
+
+
+int main()
+{
+  std::ofstream logfile("output");
+  deallog << std::setprecision(2);
+  deallog.attach(logfile);
+  deallog.depth_console(0);
+  deallog.threshold_double(1e-50); // exact equality required!
+
+  check_norms<float>();
+  check_norms<double>();
+  check_norms<long double>();
+  check_norms<std::complex<double> >();
+}
+
+
diff --git a/tests/lac/vector_accumulation_02.output b/tests/lac/vector_accumulation_02.output
new file mode 100644 (file)
index 0000000..63e9012
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
+DEAL::Vector difference: 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
diff --git a/tests/lac/vector_add_and_dot.cc b/tests/lac/vector_add_and_dot.cc
new file mode 100644 (file)
index 0000000..8b871cc
--- /dev/null
@@ -0,0 +1,79 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2012 - 2013 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 Vector::add_and_dot works correctly
+
+#include "../tests.h"
+#include <deal.II/base/logstream.h>
+#include <deal.II/lac/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;
+      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:   ";
+          v1.print(deallog);
+          deallog << "Vector check reference: ";
+          check.print(deallog);
+        }
+
+      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/vector_add_and_dot.output b/tests/lac/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

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.