From: kronbichler Date: Sun, 18 Nov 2012 08:36:09 +0000 (+0000) Subject: Parallelize inner products, norms, etc in dealii::Vector with TBB. Tests for accumula... X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e33c08d285088f76257268fc991ea7350e1719bc;p=dealii-svn.git Parallelize inner products, norms, etc in dealii::Vector with TBB. Tests for accumulation. Slight cleanup of parallel operations in Vector. git-svn-id: https://svn.dealii.org/trunk@27561 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/include/deal.II/lac/vector.h b/deal.II/include/deal.II/lac/vector.h index 5e232a756a..cd000ea746 100644 --- a/deal.II/include/deal.II/lac/vector.h +++ b/deal.II/include/deal.II/lac/vector.h @@ -16,9 +16,7 @@ #include #include #include -#include #include -#include #include #include @@ -411,7 +409,7 @@ class Vector : public Subscriptor * Copy the given vector. Resize the * present vector if necessary. */ - Vector & operator= (const Vector &c); + Vector & operator= (const Vector &v); /** * Copy the given vector. Resize the @@ -1036,17 +1034,19 @@ class Vector : public Subscriptor */ Number *val; - /* + /** * Make all other vector types * friends. */ template friend class Vector; - /* + + /** * LAPACK matrices need access to * the data. */ friend class LAPACKFullMatrix; - /* + + /** * VectorView will access the * pointer. */ @@ -1140,110 +1140,25 @@ void Vector::reinit (const unsigned int n, const bool fast) +// declare function that is implemented in vector.templates.h namespace internal { namespace Vector { - template - void set_subrange (const T s, - const unsigned int begin, - const unsigned int end, - dealii::Vector &dst) - { - if (s == T()) - std::memset ((dst.begin()+begin),0,(end-begin)*sizeof(T)); - else - std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s); - } - - template - void copy_subrange (const dealii::Vector&src, - const unsigned int begin, - const unsigned int end, - dealii::Vector &dst) - { - memcpy(&*(dst.begin()+begin), &*(src.begin()+begin), - (end-begin)*sizeof(T)); - } - - template - void copy_subrange_ext (const dealii::Vector&src, - const unsigned int begin, - const unsigned int end, - dealii::Vector &dst) - { - const T* q = src.begin()+begin; - const T* const end_q = src.begin()+end; - U* p = dst.begin()+begin; - for (; q!=end_q; ++q, ++p) - *p = *q; - } + template + void copy_vector (const dealii::Vector &src, + dealii::Vector &dst); } } -template -inline -Vector & Vector::operator = (const Number s) -{ - Assert (numbers::is_finite(s), ExcNumberNotFinite()); - if (s != Number()) - Assert (vec_size!=0, ExcEmptyObject()); - if (vec_size>dealii::internal::Vector::minimum_parallel_grain_size) - parallel::apply_to_subranges (0U, vec_size, - std_cxx1x::bind(&dealii::internal::Vector::template - set_subrange, - s, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::ref(*this)), - dealii::internal::Vector::minimum_parallel_grain_size); - else if (vec_size > 0) - dealii::internal::Vector::set_subrange(s, 0U, vec_size, *this); - - return *this; -} - - - -#ifdef DEAL_II_BOOST_BIND_COMPILER_BUG -template <> -inline -Vector > & -Vector >::operator = (const std::complex s) -{ - Assert (numbers::is_finite(s), ExcNumberNotFinite()); - if (s != std::complex()) - Assert (vec_size!=0, ExcEmptyObject()); - if (vec_size!=0) - std::fill (begin(), end(), s); - - return *this; -} -#endif - - - template inline Vector & Vector::operator = (const Vector& v) { - // if v is the same vector as *this, there is - // nothing to - if (PointerComparison::equal(this, &v) == true) - return *this; - - if (v.vec_size != vec_size) - reinit (v.vec_size, true); - if (vec_size>dealii::internal::Vector::minimum_parallel_grain_size) - parallel::apply_to_subranges (0U, vec_size, - std_cxx1x::bind(&dealii::internal::Vector::template - copy_subrange, - std_cxx1x::cref(v), std_cxx1x::_1, std_cxx1x::_2, - std_cxx1x::ref(*this)), - dealii::internal::Vector::minimum_parallel_grain_size); - else if (vec_size > 0) - dealii::internal::Vector::copy_subrange(v, 0U, vec_size, *this); - + internal::Vector::copy_vector (v, *this); return *this; } @@ -1255,18 +1170,7 @@ inline Vector & Vector::operator = (const Vector& v) { - if (v.vec_size != vec_size) - reinit (v.vec_size, true); - if (vec_size>dealii::internal::Vector::minimum_parallel_grain_size) - parallel::apply_to_subranges (0U, vec_size, - std_cxx1x::bind(&dealii::internal::Vector::template - copy_subrange_ext, - std_cxx1x::cref(v), std_cxx1x::_1, std_cxx1x::_2, - std_cxx1x::ref(*this)), - dealii::internal::Vector::minimum_parallel_grain_size); - else if (vec_size > 0) - dealii::internal::Vector::copy_subrange_ext(v, 0U, vec_size, *this); - + internal::Vector::copy_vector (v, *this); return *this; } @@ -1386,24 +1290,6 @@ Vector::operator /= (const Number factor) -template -inline -void -Vector::scale (const Number factor) -{ - Assert (numbers::is_finite(factor),ExcNumberNotFinite()); - - Assert (vec_size!=0, ExcEmptyObject()); - - parallel::transform (val, - val+vec_size, - val, - (factor*boost::lambda::_1), - dealii::internal::Vector::minimum_parallel_grain_size); -} - - - template template inline @@ -1452,51 +1338,6 @@ Vector::add (const unsigned int n_indices, -template -inline -void -Vector::add (const Number a, - const Vector& v) -{ - Assert (numbers::is_finite(a),ExcNumberNotFinite()); - - Assert (vec_size!=0, ExcEmptyObject()); - Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); - - parallel::transform (val, - val+vec_size, - v.val, - val, - (boost::lambda::_1 + a*boost::lambda::_2), - dealii::internal::Vector::minimum_parallel_grain_size); -} - - - -template -inline -void -Vector::sadd (const Number x, - const Number a, - const Vector& v) -{ - Assert (numbers::is_finite(x),ExcNumberNotFinite()); - Assert (numbers::is_finite(a),ExcNumberNotFinite()); - - Assert (vec_size!=0, ExcEmptyObject()); - Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); - - parallel::transform (val, - val+vec_size, - v.val, - val, - (x*boost::lambda::_1 + a*boost::lambda::_2), - dealii::internal::Vector::minimum_parallel_grain_size); -} - - - - template template inline diff --git a/deal.II/include/deal.II/lac/vector.templates.h b/deal.II/include/deal.II/lac/vector.templates.h index fcd4991b62..81ab9aa302 100644 --- a/deal.II/include/deal.II/lac/vector.templates.h +++ b/deal.II/include/deal.II/lac/vector.templates.h @@ -15,6 +15,8 @@ #include #include +#include +#include #include #include @@ -27,6 +29,7 @@ # include #endif +#include #include #include #include @@ -35,7 +38,6 @@ DEAL_II_NAMESPACE_OPEN -#define BLOCK_LEVEL 6 namespace internal { @@ -124,7 +126,7 @@ Vector::Vector (const Vector& v) { val = new Number[max_vec_size]; Assert (val != 0, ExcOutOfMemory()); - std::copy (v.begin(), v.end(), begin()); + *this = v; } } @@ -316,144 +318,412 @@ Vector::is_non_negative () const -template -template -Number Vector::operator * (const Vector& v) const +namespace internal { - Assert (vec_size!=0, ExcEmptyObject()); - - if (PointerComparison::equal (this, &v)) - return norm_sqr(); - - Assert (vec_size == v.size(), - ExcDimensionMismatch(vec_size, v.size())); + namespace Vector + { + template + void set_subrange (const T s, + const unsigned int begin, + const unsigned int end, + dealii::Vector &dst) + { + if (s == T()) + std::memset ((dst.begin()+begin),0,(end-begin)*sizeof(T)); + else + std::fill (&*(dst.begin()+begin), &*(dst.begin()+end), s); + } - const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), - *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), - *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); - const Number2 * Y = v.val; + template + void copy_subrange (const unsigned int begin, + const unsigned int end, + const dealii::Vector&src, + dealii::Vector &dst) + { + memcpy(&*(dst.begin()+begin), &*(src.begin()+begin), + (end-begin)*sizeof(T)); + } - // multiply the two vectors. we have to - // convert the elements of u to the type of - // the result vector. this is necessary - // because - // operator*(complex,complex) - // is not defined by default. do the - // operations block-wise with post-update. - // use three nested loops, which should make - // the roundoff error very small up to - // about 20m entries. in the - // end do extra loops with the remainders. - // this blocked algorithm has been proposed - // by Castaldo, Whaley and Chronopoulos - // (SIAM J. Sci. Comput. 31, 1156-1174, - // 2008) - while (X != X_end1) + template + void copy_subrange (const unsigned int begin, + const unsigned int end, + const dealii::Vector&src, + dealii::Vector &dst) { - sum1 = Number(); - for (unsigned int j=0; j::conjugate(*Y++)); - sum2 += sum3; - } - sum1 += sum2; - } - sum += sum1; + const T *q = src.begin()+begin; + const T *const end_q = src.begin()+end; + U *p = dst.begin()+begin; + for (; q!=end_q; ++q, ++p) + *p = *q; } - while (X != X_end2) + + template + void copy_subrange_wrap (const unsigned int begin, + const unsigned int end, + const dealii::Vector&src, + dealii::Vector &dst) { - sum2 = Number(); - for (unsigned int j=0; j::conjugate(*Y++)); - sum2 += sum3; - } - sum += sum2; + copy_subrange (begin, end, src, dst); } - while (X != X_end3) + + template + void copy_vector (const dealii::Vector&src, + dealii::Vector &dst) { - sum3 = Number(); - for (unsigned int i=0; i::conjugate(*Y++)); - sum += sum3; + const unsigned int vec_size = src.size(); + const unsigned int dst_size = dst.size(); + if (dst_size != vec_size) + dst.reinit (vec_size, true); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::apply_to_subranges (0U, vec_size, + std_cxx1x::bind(&internal::Vector::template + copy_subrange_wrap , + std_cxx1x::_1, + std_cxx1x::_2, + std_cxx1x::cref(src), + std_cxx1x::ref(dst)), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + copy_subrange (0U, vec_size, src, dst); } + } +} - sum3 = Number(); - while (X != X_end) - sum3 += *X++ * Number(numbers::NumberTraits::conjugate(*Y++)); - sum += sum3; - Assert(numbers::is_finite(sum), ExcNumberNotFinite()); - return sum; +template +Vector & +Vector::operator = (const Number s) +{ + Assert (numbers::is_finite(s), ExcNumberNotFinite()); + if (s != Number()) + Assert (vec_size!=0, ExcEmptyObject()); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::apply_to_subranges (0U, vec_size, + std_cxx1x::bind(&internal::Vector::template + set_subrange, + s, std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::ref(*this)), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + internal::Vector::set_subrange(s, 0U, vec_size, *this); + + return *this; +} + + + +#ifdef DEAL_II_BOOST_BIND_COMPILER_BUG +template <> +Vector > & +Vector >::operator = (const std::complex s) +{ + Assert (numbers::is_finite(s), ExcNumberNotFinite()); + if (s != std::complex()) + Assert (vec_size!=0, ExcEmptyObject()); + if (vec_size!=0) + std::fill (begin(), end(), s); + + return *this; } +#endif + template -typename Vector::real_type -Vector::norm_sqr () const +void +Vector::scale (const Number factor) +{ + Assert (numbers::is_finite(factor),ExcNumberNotFinite()); + + Assert (vec_size!=0, ExcEmptyObject()); + + parallel::transform (val, + val+vec_size, + val, + (factor*boost::lambda::_1), + internal::Vector::minimum_parallel_grain_size); +} + + + +template +void +Vector::add (const Number a, + const Vector& v) +{ + Assert (numbers::is_finite(a),ExcNumberNotFinite()); + + Assert (vec_size!=0, ExcEmptyObject()); + Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); + + parallel::transform (val, + val+vec_size, + v.val, + val, + (boost::lambda::_1 + a*boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); +} + + + +template +void +Vector::sadd (const Number x, + const Number a, + const Vector& v) { + Assert (numbers::is_finite(x),ExcNumberNotFinite()); + Assert (numbers::is_finite(a),ExcNumberNotFinite()); + Assert (vec_size!=0, ExcEmptyObject()); + Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); + + parallel::transform (val, + val+vec_size, + v.val, + val, + (x*boost::lambda::_1 + a*boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); +} + - const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), - *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), - *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); - while (X != X_end1) +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 + template + struct InnerProd + { + Number + operator() (const Number*&X, const Number2*&Y, const Number &) const + { + return *X++ * Number(numbers::NumberTraits::conjugate(*Y++)); + } + }; + + template + struct Norm2 + { + RealType + operator() (const Number*&X, const Number* &, const RealType &) const + { + return numbers::NumberTraits::abs_square(*X++); + } + }; + + template + struct Norm1 + { + RealType + operator() (const Number*&X, const Number* &, const RealType &) const + { + return numbers::NumberTraits::abs(*X++); + } + }; + + template + struct NormP { - sum1 = 0.; - for (unsigned int j=0; j::abs(*X++), p); + } + }; + + template + struct MeanValue + { + Number + operator() (const Number*&X, const Number* &, const Number &) const + { + return *X++; + } + }; + + // 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 + // a similar form by Castaldo, Whaley and Chronopoulos (SIAM + // J. Sci. Comput. 31, 1156-1174, 2008) and we use the smallest possible + // block size, 2. Sometimes it is referred to as pairwise summation. The + // worst case error made by this algorithm is on the order O(eps * + // log2(vec_size)), whereas a naive summation is O(eps * vec_size). Even + // though the Kahan summation is even more accurate with an error O(eps) + // by carrying along remainders not captured by the sum, that involves + // additional costs which are not worthwhile. See the Wikipedia article on + // the Kahan summation algorithm. + + // The algorithm implemented here has the additional benefit that it is + // easily parallelized without changing the order of how the elements are + // added (floating point addition is not associative). For the same vector + // size and minimum_parallel_grainsize, the blocks are always the + // same. Only at the innermost level, eight values are summed up + // consecutively in order to better balance multiplications and additions. + + // The code returns the result as the last argument in order to make + // spawning tasks simpler and use automatic template deduction. + template + void accumulate (const Operation &op, + const Number *X, + const Number2 *Y, + const ResultType power, + const unsigned int vec_size, + ResultType &result) + { + if (vec_size <= 4096) { - sum2 = 0.; - for (unsigned int k=0; k::abs_square(*X++); - sum2 += sum3; + ResultType r0 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r0 += op(X, Y, power); + ResultType r1 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r1 += op(X, Y, power); + r0 += r1; + r1 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r1 += op(X, Y, power); + ResultType r2 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r2 += op(X, Y, power); + r1 += r2; + r0 += r1; + outer_results[i] = r0; } - sum1 += sum2; + + // now work on the remainder, i.e., the last + // 32 values. Use switch statement with + // fall-through to work on these values. + if (remainder > 0) + { + const unsigned int inner_chunks = remainder / 8; + Assert (inner_chunks <= 3, ExcInternalError()); + const unsigned int remainder_inner = remainder % 8; + ResultType r0 = ResultType(), r1 = ResultType(), + r2 = ResultType(); + switch (inner_chunks) + { + case 3: + r2 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r2 += op(X, Y, power); + // no break + case 2: + r1 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r1 += op(X, Y, power); + r1 += r2; + // no break + case 1: + r2 = op(X, Y, power); + for (unsigned int j=1; j<8; ++j) + r2 += op(X, Y, power); + // no break + default: + for (unsigned int j=0; j 1) + { + if (n_chunks % 2 == 1) + outer_results[n_chunks++] = ResultType(); + for (unsigned int i=0; i 4 * internal::Vector::minimum_parallel_grain_size) { - sum3 = 0.; - for (unsigned int i=0; i::abs_square(*X++); - sum2 += sum3; + // split the vector into smaller pieces to be + // worked on recursively and create tasks for + // them. Make pieces divisible by 1024. + const unsigned int new_size = (vec_size / 4096) * 1024; + ResultType r0, r1, r2, r3; + Threads::TaskGroup<> task_group; + task_group += Threads::new_task(&accumulate, + op, X, Y, power, new_size, r0); + task_group += Threads::new_task(&accumulate, + op, X+new_size, Y+new_size, power, + new_size, r1); + task_group += Threads::new_task(&accumulate, + op, X+2*new_size, Y+2*new_size, power, + new_size, r2); + task_group += Threads::new_task(&accumulate, + op, X+3*new_size, Y+3*new_size, power, + vec_size-3*new_size, r3); + task_group.join_all(); + r0 += r1; + r2 += r3; + result = r0 + r2; + } +#endif + else + { + // split vector into four pieces and work on + // the pieces recursively. Make pieces + // divisible by 1024. + const unsigned int 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); + r0 += r1; + r2 += r3; + result = r0 + r2; } - sum += sum2; - } - while (X != X_end3) - { - sum3 = 0.; - for (unsigned int i=0; i::abs_square(*X++); - sum += sum3; } + } +} - sum3 = 0.; - while (X != X_end) - sum3 += numbers::NumberTraits::abs_square(*X++); - sum += sum3; + +template +template +Number Vector::operator * (const Vector& v) const +{ + Assert (vec_size!=0, ExcEmptyObject()); + + if (PointerComparison::equal (this, &v)) + return norm_sqr(); + + Assert (vec_size == v.size(), + ExcDimensionMismatch(vec_size, v.size())); + + Number sum; + internal::Vector::accumulate (internal::Vector::InnerProd(), + val, v.val, Number(), vec_size, sum); Assert(numbers::is_finite(sum), ExcNumberNotFinite()); return sum; @@ -461,58 +731,29 @@ Vector::norm_sqr () const template -Number Vector::mean_value () const +typename Vector::real_type +Vector::norm_sqr () const { Assert (vec_size!=0, ExcEmptyObject()); - const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), - *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), - *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); + real_type sum; + internal::Vector::accumulate (internal::Vector::Norm2(), + val, val, real_type(), vec_size, sum); - while (X != X_end1) - { - sum1 = Number(); - for (unsigned int j=0; j +Number Vector::mean_value () const +{ + Assert (vec_size!=0, ExcEmptyObject()); + + Number sum; + internal::Vector::accumulate (internal::Vector::MeanValue(), + val, val, Number(), vec_size, sum); return sum / real_type(size()); } @@ -525,54 +766,9 @@ Vector::l1_norm () const { Assert (vec_size!=0, ExcEmptyObject()); - const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), - *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), - *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); - - while (X != X_end1) - { - sum1 = 0.; - for (unsigned int j=0; j::abs(*X++); - sum2 += sum3; - } - sum1 += sum2; - } - sum += sum1; - } - while (X != X_end2) - { - sum2 = 0.; - for (unsigned int j=0; j::abs(*X++); - sum2 += sum3; - } - sum += sum2; - } - while (X != X_end3) - { - sum3 = 0.; - for (unsigned int i=0; i::abs(*X++); - sum += sum3; - } - - sum3 = 0.; - while (X != X_end) - sum3 += numbers::NumberTraits::abs(*X++); - sum += sum3; + real_type sum; + internal::Vector::accumulate (internal::Vector::Norm1(), + val, val, real_type(), vec_size, sum); return sum; } @@ -583,7 +779,41 @@ template typename Vector::real_type Vector::l2_norm () const { - return std::sqrt(norm_sqr()); + // if l2_norm()^2 is finite and non-zero, the answer is computed as + // std::sqrt(norm_sqr()). If norm_sqr() is infinite or zero, the l2 norm + // might still be finite. In that case, recompute it (this is a rare case, + // so working on the vector twice is uncritical and paid off by the extended + // precision) using the BLAS approach with a weight, see e.g. dnrm2.f. + Assert (vec_size!=0, ExcEmptyObject()); + + real_type norm_square; + internal::Vector::accumulate (internal::Vector::Norm2(), + val, val, real_type(), vec_size, norm_square); + if (numbers::is_finite(norm_square) && + norm_square > std::numeric_limits::min()) + return std::sqrt(norm_square); + else + { + real_type scale = 0.; + real_type sum = 1.; + for (unsigned int i=0; i::abs(val[i]); + if (scale < abs_x) + { + sum = 1. + sum * (scale/abs_x) * (scale/abs_x); + scale = abs_x; + } + else + sum += (abs_x/scale) * (abs_x/scale); + } + } + Assert(numbers::is_finite(scale)*std::sqrt(sum), ExcNumberNotFinite()); + return scale * std::sqrt(sum); + } } @@ -597,58 +827,35 @@ Vector::lp_norm (const real_type p) const if (p == 1.) return l1_norm(); else if (p == 2.) - return std::sqrt(norm_sqr()); + return l2_norm(); - const unsigned int blocking = 1<>(BLOCK_LEVEL))<<(BLOCK_LEVEL)), - *X_end2 = X + ((vec_size>>(2*BLOCK_LEVEL))<<(2*BLOCK_LEVEL)), - *X_end1 = X + ((vec_size>>(3*BLOCK_LEVEL))<<(3*BLOCK_LEVEL)); + real_type sum; + internal::Vector::accumulate (internal::Vector::NormP(), + val, val, p, vec_size, sum); - while (X != X_end1) + if (numbers::is_finite(sum) && sum > std::numeric_limits::min()) + return std::pow(sum, static_cast(1./p)); + else { - sum1 = 0.; - for (unsigned int j=0; j::abs(*X++),p); - sum2 += sum3; + const real_type abs_x = + numbers::NumberTraits::abs(val[i]); + if (scale < abs_x) + { + sum = 1. + sum * std::pow(scale/abs_x, p); + scale = abs_x; + } + else + sum += std::pow(abs_x/scale, p); } - sum1 += sum2; - } - sum += sum1; - } - while (X != X_end2) - { - sum2 = 0.; - for (unsigned int j=0; j::abs(*X++),p); - sum2 += sum3; } - sum += sum2; + return scale * std::pow(sum, static_cast(1./p)); } - while (X != X_end3) - { - sum3 = 0.; - for (unsigned int i=0; i::abs(*X++),p); - sum += sum3; - } - - sum3 = 0.; - while (X != X_end) - sum3 += std::pow(numbers::NumberTraits::abs(*X++),p); - sum += sum3; - - return std::pow(sum, static_cast(1./p)); } @@ -797,13 +1004,13 @@ void Vector::sadd (const Number x, const Number a, Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size)); if (vec_size>internal::Vector::minimum_parallel_grain_size) - parallel::transform (val, - val+vec_size, - v.val, - w.val, - val, - (x*boost::lambda::_1 + a*boost::lambda::_2 + b*boost::lambda::_3), - internal::Vector::minimum_parallel_grain_size); + parallel::transform (val, + val+vec_size, + v.val, + w.val, + val, + (x*boost::lambda::_1 + a*boost::lambda::_2 + b*boost::lambda::_3), + internal::Vector::minimum_parallel_grain_size); else if (vec_size > 0) for (unsigned int i=0; i::operator = (const PETScWrappers::Vector &v) internal::copy (start_ptr, start_ptr+vec_size, begin()); - // restore the representation of the - // vector - ierr = VecRestoreArray (static_cast(v), &start_ptr); + // restore the representation of the + // vector + ierr = VecRestoreArray (static_cast(v), &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); } @@ -1145,7 +1352,7 @@ void Vector::print (std::ostream &out, template void -Vector::print (LogStream& out, const unsigned int width, const bool across) const +Vector::print (LogStream &out, const unsigned int width, const bool across) const { Assert (vec_size!=0, ExcEmptyObject()); @@ -1177,11 +1384,11 @@ void Vector::block_write (std::ostream &out) const std::strcat(buf, "\n["); out.write(buf, std::strlen(buf)); - out.write (reinterpret_cast(begin()), - reinterpret_cast(end()) - - reinterpret_cast(begin())); + out.write (reinterpret_cast(begin()), + reinterpret_cast(end()) + - reinterpret_cast(begin())); - // out << ']'; + // out << ']'; const char outro = ']'; out.write (&outro, 1); @@ -1208,15 +1415,15 @@ void Vector::block_read (std::istream &in) reinit (sz, true); char c; - // in >> c; + // in >> c; in.read (&c, 1); AssertThrow (c=='[', ExcIO()); - in.read (reinterpret_cast(begin()), - reinterpret_cast(end()) - - reinterpret_cast(begin())); + in.read (reinterpret_cast(begin()), + reinterpret_cast(end()) + - reinterpret_cast(begin())); - // in >> c; + // in >> c; in.read (&c, 1); AssertThrow (c==']', ExcIO()); } diff --git a/deal.II/source/grid/grid_reordering.cc b/deal.II/source/grid/grid_reordering.cc index 429490e4f6..3355f06ded 100644 --- a/deal.II/source/grid/grid_reordering.cc +++ b/deal.II/source/grid/grid_reordering.cc @@ -15,6 +15,7 @@ #include #include #include +#include #include #include diff --git a/deal.II/source/grid/tria.cc b/deal.II/source/grid/tria.cc index ae9148a0e0..eb5dc1836f 100644 --- a/deal.II/source/grid/tria.cc +++ b/deal.II/source/grid/tria.cc @@ -13,6 +13,8 @@ #include #include +#include +#include #include #include @@ -20,7 +22,6 @@ #include #include #include -#include #include #include #include diff --git a/deal.II/source/hp/dof_handler.cc b/deal.II/source/hp/dof_handler.cc index 1e162dcfb9..ce43f0da7a 100644 --- a/deal.II/source/hp/dof_handler.cc +++ b/deal.II/source/hp/dof_handler.cc @@ -13,6 +13,7 @@ #include +#include #include #include #include diff --git a/deal.II/source/lac/vector.inst.in b/deal.II/source/lac/vector.inst.in index 2ac1a71afb..73bdab2d77 100644 --- a/deal.II/source/lac/vector.inst.in +++ b/deal.II/source/lac/vector.inst.in @@ -19,6 +19,14 @@ for (SCALAR : REAL_SCALARS) for (S1, S2 : REAL_SCALARS) { + namespace internal + \{ + namespace Vector + \{ + template void copy_vector (const dealii::Vector&, + dealii::Vector&); + \} + \} template bool Vector::operator==(const Vector&) const; @@ -38,6 +46,14 @@ for (SCALAR : COMPLEX_SCALARS) for (S1, S2 : COMPLEX_SCALARS) { + namespace internal + \{ + namespace Vector + \{ + template void copy_vector (const dealii::Vector&, + dealii::Vector&); + \} + \} template bool Vector::operator==(const Vector&) const; diff --git a/tests/lac/vector_accumulation.cc b/tests/lac/vector_accumulation.cc new file mode 100644 index 0000000000..d4f6bf8294 --- /dev/null +++ b/tests/lac/vector_accumulation.cc @@ -0,0 +1,64 @@ +//----------------------- vector_accumulation.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//----------------------- vector_accumulation.cc --------------------------- + +// check that the l2 norm is exactly the same for many runs on random vector +// size with random vector entries. this is to ensure that the result is +// reproducible also when parallel evaluation is done + +#include "../tests.h" +#include +#include +#include +#include +#include + + + + +template +void check_norms () +{ + for (unsigned int test=0; test<20; ++test) + { + const unsigned int size = rand() % 100000; + Vector vec (size); + for (unsigned int i=0; i(rand())/static_cast(RAND_MAX); + const typename Vector::real_type norm = vec.l2_norm(); + for (unsigned int i=0; i<30; ++i) + Assert (vec.l2_norm() == norm, ExcInternalError()); + + Vector vec2 (vec); + for (unsigned int i=0; i<10; ++i) + Assert (vec2.l2_norm() == norm, ExcInternalError()); + } +} + + +int main() +{ + std::ofstream logfile("vector_accumulation/output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check_norms(); + check_norms(); + check_norms(); + check_norms >(); + deallog << "OK" << std::endl; +} + + diff --git a/tests/lac/vector_accumulation/cmp/generic b/tests/lac/vector_accumulation/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/vector_accumulation/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/vector_large_numbers.cc b/tests/lac/vector_large_numbers.cc new file mode 100644 index 0000000000..edb5b8f4f3 --- /dev/null +++ b/tests/lac/vector_large_numbers.cc @@ -0,0 +1,77 @@ +//------------------------ vector_large_numbers.cc ------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//------------------------ vector_large_numbers.cc ------------------------- + + +#include "../tests.h" +#include +#include +#include +#include +#include + + + +void check_large_numbers() +{ + Vector v(10); + v(0) = 1e13; + v(1) = 1e19; + v(2) = 1e21; + v(4) = 1.; + v(7) = 1e20; + + const double correct = std::sqrt(1e26 + 1e38 + 1e42 + 1e40 + 1.); + Assert (std::abs(v.l2_norm() - correct) < 1e-6*correct, ExcInternalError()); + + v = 0.; + v(5) = 1e-30; + v(9) = 1e-32; + const double correct2 = std::sqrt(1e-60 + 1e-64); + Assert (std::abs(v.l2_norm() - correct2) < 1e-6*correct2, ExcInternalError()); + + Vector w(7); + w(0) = 1e232; + w(6) = 1e231; + w(3) = 1e200; + const double correct3 = std::sqrt(100. + 1.) * 1e231; + Assert (std::abs(w.l2_norm() - correct3) < 1e-13*correct3, ExcInternalError()); + + w = 0; + w(1) = 1e-302; + w(2) = 1e-303; + w(3) = 2e-303; + w(4) = 3e-303; + w(5) = -1e-303; + const double correct4 = std::sqrt(100. + 1. + 4. + 9 + 1.) * 1e-303; + Assert (std::abs(w.l2_norm() - correct4) < 1e-13*correct4, ExcInternalError()); + + const double correct5 = std::pow(1000. + 1. + 8. + 27 + 1., 1./3.) * 1e-303; + Assert (std::abs(w.lp_norm(3.) - correct5) < 1e-13*correct5, ExcInternalError()); + deallog << "OK" << std::endl; +} + + + +int main() +{ + std::ofstream logfile("vector_large_numbers/output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check_large_numbers(); +} + + diff --git a/tests/lac/vector_large_numbers/cmp/generic b/tests/lac/vector_large_numbers/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/vector_large_numbers/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK diff --git a/tests/lac/vector_norms.cc b/tests/lac/vector_norms.cc new file mode 100644 index 0000000000..9ada6fa924 --- /dev/null +++ b/tests/lac/vector_norms.cc @@ -0,0 +1,126 @@ +//---------------------------- vector_norms.cc --------------------------- +// $Id$ +// Version: $Name$ +// +// Copyright (C) 2012 by the deal.II authors +// +// This file is subject to QPL and may not be distributed +// without copyright and license information. Please refer +// to the file deal.II/doc/license.html for the text and +// further information on this license. +// +//---------------------------- vector_norms.cc --------------------------- + +// checks that the l1, l2, lp norm is computed correctly for deal.II vectors +// (check the summation algorithm), including an accuracy test (should not +// lose more than 1 decimal also for 200000 vector entries) + +#include "../tests.h" +#include +#include +#include +#include +#include + + + + +template +void check_norms () +{ + const number acc = 1e1*std::numeric_limits::epsilon(); + unsigned int skip = 73; + for (unsigned int size=1; size<200000; size+=skip) + { + // test correctness + if (size > 10000) + skip += 17; + Vector vec(size); + for (unsigned int i=0; i(size), + static_cast(1./3.))) < + std::max(acc,number(1e-15))*std::pow(static_cast(size), + static_cast(1./3.))*value, + ExcInternalError()); + } +} + + +template +void check_complex_norms () +{ + const number acc = 1e2*std::numeric_limits::epsilon(); + unsigned int skip = 73; + for (unsigned int size=1; size<100000; size+=skip) + { + // test correctness + if (size > 10000) + skip += 17; + Vector > vec(size); + long double sum = 0.; + for (unsigned int i=0; i(i+1,i+2); + sum += std::sqrt((long double)(i+1)*(i+1) + (long double)(i+2)*(i+2)); + } + + const number l1_norm = vec.l1_norm(); + Assert (std::abs(l1_norm-sum) < acc*sum, + ExcInternalError()); + + // test accuracy of summation + const std::complex value (3.14159265358979323846, 0.1); + vec = std::complex(value); + const number l1_norma = vec.l1_norm(); + Assert (std::abs(l1_norma-std::abs(value)*size) < + acc*size*std::abs(value), + ExcInternalError()); + const number l2_norma = vec.l2_norm(); + Assert (std::abs(l2_norma-std::abs(value)*std::sqrt((number)size)) < + acc*std::sqrt((number)size)*std::abs(value), + ExcInternalError()); + const number lp_norma = vec.lp_norm(3.); + Assert (std::abs(lp_norma-std::abs(value)* + std::pow(static_cast(size), + static_cast(1./3.))) < + acc*std::pow(static_cast(size), + static_cast(1./3.))*std::abs(value), + ExcInternalError()); + } +} + + +int main() +{ + std::ofstream logfile("vector_norms/output"); + deallog << std::fixed; + deallog << std::setprecision(2); + deallog.attach(logfile); + deallog.depth_console(0); + deallog.threshold_double(1.e-10); + + check_norms(); + check_norms(); + check_norms(); + check_complex_norms(); + check_complex_norms(); + deallog << "OK" << std::endl; +} + + diff --git a/tests/lac/vector_norms/cmp/generic b/tests/lac/vector_norms/cmp/generic new file mode 100644 index 0000000000..0fd8fc12f0 --- /dev/null +++ b/tests/lac/vector_norms/cmp/generic @@ -0,0 +1,2 @@ + +DEAL::OK