From abcd2880e668632c64ab6797ecbba52ff994e78a Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Sun, 26 Oct 2014 09:25:43 +0100 Subject: [PATCH] Combined add and inner product for vectors: add_and_dot - 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 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. --- doc/news/changes.h | 15 + include/deal.II/base/vectorization.h | 23 ++ include/deal.II/lac/block_vector_base.h | 40 +++ include/deal.II/lac/parallel_block_vector.h | 43 ++- include/deal.II/lac/parallel_vector.h | 65 +++- .../deal.II/lac/parallel_vector.templates.h | 10 +- include/deal.II/lac/petsc_vector_base.h | 19 ++ include/deal.II/lac/solver_bicgstab.h | 12 +- include/deal.II/lac/solver_cg.h | 3 +- include/deal.II/lac/solver_gmres.h | 41 +-- include/deal.II/lac/trilinos_vector_base.h | 31 ++ include/deal.II/lac/vector.h | 25 ++ include/deal.II/lac/vector.templates.h | 295 ++++++++++++++---- source/lac/petsc_vector_base.cc | 13 +- tests/lac/gmres_reorthogonalize_02.cc | 4 +- ...nalize_02.compiler=Intel=no.release.output | 5 - ...alize_02.compiler=Intel=yes.release.output | 5 - .../lac/gmres_reorthogonalize_02.debug.output | 5 - tests/lac/gmres_reorthogonalize_02.output | 5 + tests/lac/gmres_reorthogonalize_04.cc | 2 +- ...nalize_04.compiler=Intel=no.release.output | 5 - .../lac/gmres_reorthogonalize_04.debug.output | 5 - ...output => gmres_reorthogonalize_04.output} | 2 +- ...umulation.cc => vector_accumulation_01.cc} | 0 ...n.output => vector_accumulation_01.output} | 0 tests/lac/vector_accumulation_02.cc | 82 +++++ tests/lac/vector_accumulation_02.output | 21 ++ tests/lac/vector_add_and_dot.cc | 79 +++++ tests/lac/vector_add_and_dot.output | 24 ++ 29 files changed, 759 insertions(+), 120 deletions(-) delete mode 100644 tests/lac/gmres_reorthogonalize_02.compiler=Intel=no.release.output delete mode 100644 tests/lac/gmres_reorthogonalize_02.compiler=Intel=yes.release.output delete mode 100644 tests/lac/gmres_reorthogonalize_02.debug.output create mode 100644 tests/lac/gmres_reorthogonalize_02.output delete mode 100644 tests/lac/gmres_reorthogonalize_04.compiler=Intel=no.release.output delete mode 100644 tests/lac/gmres_reorthogonalize_04.debug.output rename tests/lac/{gmres_reorthogonalize_04.compiler=Intel=yes.release.output => gmres_reorthogonalize_04.output} (70%) rename tests/lac/{vector_accumulation.cc => vector_accumulation_01.cc} (100%) rename tests/lac/{vector_accumulation.output => vector_accumulation_01.output} (100%) create mode 100644 tests/lac/vector_accumulation_02.cc create mode 100644 tests/lac/vector_accumulation_02.output create mode 100644 tests/lac/vector_add_and_dot.cc create mode 100644 tests/lac/vector_add_and_dot.output diff --git a/doc/news/changes.h b/doc/news/changes.h index ce74a29e2b..eb0567c0d3 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -307,6 +307,21 @@ inconvenience this causes.

Specific improvements

    +
  1. 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. +
    + (Martin Kronbichler, 2014/10/27) +
    2. + +
  2. Improved: Inner products and norms on deal.II's own vector classes now + use vectorization through VectorizedArray if available. +
    + (Martin Kronbichler, 2014/10/27) +
    3. +
  3. Changed: PETSc and Trilinos vectors with ghost entries can now be reset to zero using = 0.0;
    diff --git a/include/deal.II/base/vectorization.h b/include/deal.II/base/vectorization.h index d0993978b9..7b903e19db 100644 --- a/include/deal.II/base/vectorization.h +++ b/include/deal.II/base/vectorization.h @@ -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 + * {pow(x[0],p), pow(x[1],p), ..., pow(x[n_array_elements-1],p)}. + * + * @relates VectorizedArray + */ + template + inline + ::dealii::VectorizedArray + pow (const ::dealii::VectorizedArray &x, + const Number p) + { + Number values[::dealii::VectorizedArray::n_array_elements]; + for (unsigned int i=0; i::n_array_elements; ++i) + values[i] = std::pow(x[i], p); + ::dealii::VectorizedArray 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 {abs(x[0]), diff --git a/include/deal.II/lac/block_vector_base.h b/include/deal.II/lac/block_vector_base.h index 2a8f345808..543b4a495e 100644 --- a/include/deal.II/lac/block_vector_base.h +++ b/include/deal.II/lac/block_vector_base.h @@ -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 and parallel::distributed::Vector). + * 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::linfty_norm () const +template +typename BlockVectorBase::value_type +BlockVectorBase:: +add_and_dot (const typename BlockVectorBase::value_type a, + const BlockVectorBase &V, + const BlockVectorBase &W) +{ + AssertDimension (n_blocks(), V.n_blocks()); + AssertDimension (n_blocks(), W.n_blocks()); + + value_type sum = 0.; + for (size_type i=0; i BlockVectorBase & BlockVectorBase::operator += (const BlockVectorBase &v) diff --git a/include/deal.II/lac/parallel_block_vector.h b/include/deal.II/lac/parallel_block_vector.h index 2c38416dfb..2ca7098814 100644 --- a/include/deal.II/lac/parallel_block_vector.h +++ b/include/deal.II/lac/parallel_block_vector.h @@ -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 &V, + const BlockVector &W); + /** * Scale each element of the vector by the given factor. * @@ -928,7 +948,28 @@ namespace parallel template inline - void BlockVector::swap (BlockVector &v) + Number + BlockVector::add_and_dot (const Number a, + const BlockVector &V, + const BlockVector &W) + { + Number local_result = Number(); + for (unsigned int i=0; in_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 + inline + void + BlockVector::swap (BlockVector &v) { Assert (this->n_blocks() == v.n_blocks(), ExcDimensionMismatch(this->n_blocks(), v.n_blocks())); diff --git a/include/deal.II/lac/parallel_vector.h b/include/deal.II/lac/parallel_vector.h index 416d615f1d..06c64f01ec 100644 --- a/include/deal.II/lac/parallel_vector.h +++ b/include/deal.II/lac/parallel_vector.h @@ -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 &V, + const Vector &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 &V, + const Vector &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::~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 + inline + Number + Vector::add_and_dot_local(const Number a, + const Vector &V, + const Vector &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 + inline + Number + Vector::add_and_dot (const Number a, + const Vector &V, + const Vector &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 inline typename Vector::size_type diff --git a/include/deal.II/lac/parallel_vector.templates.h b/include/deal.II/lac/parallel_vector.templates.h index 6b2c554222..bd5b84af25 100644 --- a/include/deal.II/lac/parallel_vector.templates.h +++ b/include/deal.II/lac/parallel_vector.templates.h @@ -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 #include +#include + 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(_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; } diff --git a/include/deal.II/lac/petsc_vector_base.h b/include/deal.II/lac/petsc_vector_base.h index efb6b08201..8a67bf503e 100644 --- a/include/deal.II/lac/petsc_vector_base.h +++ b/include/deal.II/lac/petsc_vector_base.h @@ -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. diff --git a/include/deal.II/lac/solver_bicgstab.h b/include/deal.II/lac/solver_bicgstab.h index 5041c67c19..f86d60e47e 100644 --- a/include/deal.II/lac/solver_bicgstab.h +++ b/include/deal.II/lac/solver_bicgstab.h @@ -351,11 +351,11 @@ SolverBicgstab::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::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); diff --git a/include/deal.II/lac/solver_cg.h b/include/deal.II/lac/solver_cg.h index 9f57f28127..9facc394ea 100644 --- a/include/deal.II/lac/solver_cg.h +++ b/include/deal.II/lac/solver_cg.h @@ -391,9 +391,8 @@ SolverCG::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); diff --git a/include/deal.II/lac/solver_gmres.h b/include/deal.II/lac/solver_gmres.h index 9e4d933325..a5c0622ed6 100644 --- a/include/deal.II/lac/solver_gmres.h +++ b/include/deal.II/lac/solver_gmres.h @@ -528,6 +528,7 @@ SolverGMRES::modified_gram_schmidt (const internal::SolverGMRES::TmpVect Vector &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::modified_gram_schmidt (const internal::SolverGMRES::TmpVect norm_vv_start = vv.l2_norm(); // Orthogonalization - for (unsigned int i=0 ; i::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::epsilon())) return norm_vv; @@ -565,14 +564,18 @@ SolverGMRES::modified_gram_schmidt (const internal::SolverGMRES::TmpVect } if (re_orthogonalize == true) - for (unsigned int i=0 ; i::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 diff --git a/include/deal.II/lac/trilinos_vector_base.h b/include/deal.II/lac/trilinos_vector_base.h index 449c4a5cf0..5fbdf728ec 100644 --- a/include/deal.II/lac/trilinos_vector_base.h +++ b/include/deal.II/lac/trilinos_vector_base.h @@ -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 diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index 47765ee3dc..4c037a8e3b 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -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 &V, + const Vector &W); + //@} diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index d6073f5030..a331a19a83 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -22,6 +22,7 @@ #include #include #include +#include #include #include @@ -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 - 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::conjugate(*Y++)); } + + VectorizedArray + do_vectorized(const Number *&X, const Number *&Y, const Number &, Number *&) const + { + VectorizedArray x, y; + x.load(X); + y.load(Y); + X += VectorizedArray::n_array_elements; + Y += VectorizedArray::n_array_elements; + return x * y; + } }; template 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::abs_square(*X++); } + + VectorizedArray + do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const + { + VectorizedArray x; + x.load(X); + X += VectorizedArray::n_array_elements; + return x * x; + } }; template 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::abs(*X++); } + + VectorizedArray + do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const + { + VectorizedArray x; + x.load(X); + X += VectorizedArray::n_array_elements; + return std::abs(x); + } }; template 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::abs(*X++), p); } + + VectorizedArray + do_vectorized(const Number *&X, const Number *&, const Number &p, Number *&) const + { + VectorizedArray x; + x.load(X); + X += VectorizedArray::n_array_elements; + return std::pow(std::abs(x),p); + } }; template 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 + do_vectorized(const Number *&X, const Number *&, const Number &, Number *&) const + { + VectorizedArray x; + x.load(X); + X += VectorizedArray::n_array_elements; + return x; + } }; + template + struct AddAndDot + { + Number + operator() (const Number *&V, const Number *&W, const Number &a, + Number *&X) const + { + *X += a **V++; + return *X++ * Number(numbers::NumberTraits::conjugate(*W++)); + } + + VectorizedArray + do_vectorized(const Number *&V, const Number *&W, const Number &a, + Number *&X) const + { + VectorizedArray 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::n_array_elements; + V += VectorizedArray::n_array_elements; + W += VectorizedArray::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 + 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, + const unsigned int start_chunk=0) + { + AssertIndexRange(start_chunk, n_chunks+1); + for (size_type i=start_chunk; i + 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) + { + for (size_type i=0; i::n_array_elements; ++i) + { + VectorizedArray r0 = op.do_vectorized(X, Y, power, Z); + VectorizedArray r1 = op.do_vectorized(X, Y, power, Z); + VectorizedArray r2 = op.do_vectorized(X, Y, power, Z); + VectorizedArray 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::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::n_array_elements != 0) + accumulate_regular(op, X, Y, power, n_chunks, Z, outer_results, + internal::bool2type(), + (n_chunks/VectorizedArray::n_array_elements) + * VectorizedArray::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::value && + (types_are_equal::value || + types_are_equal::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 task_group; task_group += Threads::new_task(&accumulate, - 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, 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, 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, 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::operator * (const Vector &v) const ExcDimensionMismatch(vec_size, v.size())); Number sum; - internal::Vector::accumulate (internal::Vector::InnerProd(), - val, v.val, Number(), vec_size, sum); + internal::Vector::accumulate (internal::Vector::Dot(), + val, v.val, Number(), vec_size, val, sum); Assert(numbers::is_finite(sum), ExcNumberNotFinite()); return sum; } + template typename Vector::real_type Vector::norm_sqr () const @@ -1079,7 +1238,7 @@ Vector::norm_sqr () const real_type sum; internal::Vector::accumulate (internal::Vector::Norm2(), - 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::norm_sqr () const } + template Number Vector::mean_value () const { @@ -1094,7 +1254,7 @@ Number Vector::mean_value () const Number sum; internal::Vector::accumulate (internal::Vector::MeanValue(), - val, val, Number(), vec_size, sum); + val, val, Number(), vec_size, val, sum); return sum / real_type(size()); } @@ -1109,7 +1269,7 @@ Vector::l1_norm () const real_type sum; internal::Vector::accumulate (internal::Vector::Norm1(), - val, val, real_type(), vec_size, sum); + val, val, real_type(), vec_size, val, sum); return sum; } @@ -1129,7 +1289,7 @@ Vector::l2_norm () const real_type norm_square; internal::Vector::accumulate (internal::Vector::Norm2(), - 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::min()) return std::sqrt(norm_square); @@ -1172,7 +1332,7 @@ Vector::lp_norm (const real_type p) const real_type sum; internal::Vector::accumulate (internal::Vector::NormP(), - val, val, p, vec_size, sum); + val, val, p, vec_size, val, sum); if (numbers::is_finite(sum) && sum >= std::numeric_limits::min()) return std::pow(sum, static_cast(1./p)); @@ -1216,6 +1376,27 @@ Vector::linfty_norm () const } + +template +Number +Vector::add_and_dot (const Number a, + const Vector &V, + const Vector &W) +{ + Assert (vec_size!=0, ExcEmptyObject()); + AssertDimension (vec_size, V.size()); + AssertDimension (vec_size, W.size()); + + Number sum; + internal::Vector::accumulate (internal::Vector::AddAndDot(), + V.val, W.val, a, vec_size, val, sum); + Assert(numbers::is_finite(sum), ExcNumberNotFinite()); + + return sum; +} + + + template Vector &Vector::operator += (const Vector &v) { @@ -1226,6 +1407,7 @@ Vector &Vector::operator += (const Vector &v) } + template Vector &Vector::operator -= (const Vector &v) { @@ -1241,6 +1423,7 @@ Vector &Vector::operator -= (const Vector &v) } + template void Vector::add (const Number v) { @@ -1253,6 +1436,7 @@ void Vector::add (const Number v) } + template void Vector::add (const Vector &v) { @@ -1266,6 +1450,7 @@ void Vector::add (const Vector &v) } + template void Vector::add (const Number a, const Vector &v, const Number b, const Vector &w) diff --git a/source/lac/petsc_vector_base.cc b/source/lac/petsc_vector_base.cc index 041648c52e..26324e95f6 100644 --- a/source/lac/petsc_vector_base.cc +++ b/source/lac/petsc_vector_base.cc @@ -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()); diff --git a/tests/lac/gmres_reorthogonalize_02.cc b/tests/lac/gmres_reorthogonalize_02.cc index f209d1462e..bb31013603 100644 --- a/tests/lac/gmres_reorthogonalize_02.cc +++ b/tests/lac/gmres_reorthogonalize_02.cc @@ -40,7 +40,7 @@ void test () for (unsigned int i=0; i::epsilon()); + SolverControl control(1000, 1e3*std::numeric_limits::epsilon()); typename SolverGMRES >::AdditionalData data; data.max_n_tmp_vectors = 202; @@ -59,7 +59,7 @@ int main() deallog.push("double"); test(); deallog.pop(); - deallog.threshold_double(1.e-4); + deallog.threshold_double(2.e-4); deallog.push("float"); test(); 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 index 561b5e490e..0000000000 --- a/tests/lac/gmres_reorthogonalize_02.compiler=Intel=no.release.output +++ /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 index e187d5b60d..0000000000 --- a/tests/lac/gmres_reorthogonalize_02.compiler=Intel=yes.release.output +++ /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 index 561b5e490e..0000000000 --- a/tests/lac/gmres_reorthogonalize_02.debug.output +++ /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 index 0000000000..9b62b4faa2 --- /dev/null +++ b/tests/lac/gmres_reorthogonalize_02.output @@ -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 diff --git a/tests/lac/gmres_reorthogonalize_04.cc b/tests/lac/gmres_reorthogonalize_04.cc index b021e983fe..344671676a 100644 --- a/tests/lac/gmres_reorthogonalize_04.cc +++ b/tests/lac/gmres_reorthogonalize_04.cc @@ -87,7 +87,7 @@ int main() deallog.push("double"); test(); deallog.pop(); - deallog.threshold_double(1.e-4); + deallog.threshold_double(2.e-4); deallog.push("float"); test(); 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 index 06bb1801e3..0000000000 --- a/tests/lac/gmres_reorthogonalize_04.compiler=Intel=no.release.output +++ /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 index 06bb1801e3..0000000000 --- a/tests/lac/gmres_reorthogonalize_04.debug.output +++ /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.compiler=Intel=yes.release.output b/tests/lac/gmres_reorthogonalize_04.output similarity index 70% rename from tests/lac/gmres_reorthogonalize_04.compiler=Intel=yes.release.output rename to tests/lac/gmres_reorthogonalize_04.output index 84f47a4765..5d6219e851 100644 --- a/tests/lac/gmres_reorthogonalize_04.compiler=Intel=yes.release.output +++ b/tests/lac/gmres_reorthogonalize_04.output @@ -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.cc b/tests/lac/vector_accumulation_01.cc similarity index 100% rename from tests/lac/vector_accumulation.cc rename to tests/lac/vector_accumulation_01.cc diff --git a/tests/lac/vector_accumulation.output b/tests/lac/vector_accumulation_01.output similarity index 100% rename from tests/lac/vector_accumulation.output rename to tests/lac/vector_accumulation_01.output diff --git a/tests/lac/vector_accumulation_02.cc b/tests/lac/vector_accumulation_02.cc new file mode 100644 index 0000000000..9dfa75567a --- /dev/null +++ b/tests/lac/vector_accumulation_02.cc @@ -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 +#include +#include +#include +#include +#include + + + + +template +void check_norms () +{ + for (unsigned int test=0; test<5; ++test) + { + const unsigned int size = Testing::rand() % 20000; + Vector larger1 (size+8), larger2(size+8), in1(size), in2(size); + for (unsigned int i=0; i(Testing::rand())/static_cast(RAND_MAX); + in2(i) = static_cast(Testing::rand())/static_cast(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 v1 (size, larger1.begin()+shift1); + VectorView v2 (size, larger2.begin()+shift2); + for (unsigned int i=0; i(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(); + check_norms(); + check_norms(); + check_norms >(); +} + + diff --git a/tests/lac/vector_accumulation_02.output b/tests/lac/vector_accumulation_02.output new file mode 100644 index 0000000000..63e9012248 --- /dev/null +++ b/tests/lac/vector_accumulation_02.output @@ -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 index 0000000000..8b871cc8dc --- /dev/null +++ b/tests/lac/vector_add_and_dot.cc @@ -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 +#include +#include +#include +#include + + + + +template +void check () +{ + for (unsigned int test=0; test<5; ++test) + { + const unsigned int size = 17 + test*1101; + Vector v1 (size), v2(size), v3(size), check(size); + for (unsigned int i=0; i::value) + { + deallog << "Vector add reference: "; + v1.print(deallog); + deallog << "Vector check reference: "; + check.print(deallog); + } + + deallog << "Add and dot should be " << prod/static_cast(size) + << ", is " << prod_check/static_cast(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(); + check(); + check(); + check >(); + 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 index 0000000000..329e9f8d54 --- /dev/null +++ b/tests/lac/vector_add_and_dot.output @@ -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 -- 2.39.5