From: Martin Kronbichler Date: Thu, 5 Nov 2009 01:14:00 +0000 (+0000) Subject: Change to an implementation that has reduced roundoff error accumulation for the... X-Git-Tag: v8.0.0~6856 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=17d5a5b4d798d53c5a02f8dcef410f047e47d257;p=dealii.git Change to an implementation that has reduced roundoff error accumulation for the inner products and norms. Do this by using nested blocking with post-update of the sum. Also change some of the parallelizations by not going through TBB in case we work with small vectors. git-svn-id: https://svn.dealii.org/trunk@20028 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/lac/include/lac/vector.h b/deal.II/lac/include/lac/vector.h index 9cade6f7f1..275d8437df 100644 --- a/deal.II/lac/include/lac/vector.h +++ b/deal.II/lac/include/lac/vector.h @@ -1103,16 +1103,19 @@ inline Vector & Vector::operator = (const Number s) { Assert (numbers::is_finite(s), - ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); + ExcMessage("The given value is not finite but either infinite or " + "Not A Number (NaN)")); if (s != Number()) Assert (vec_size!=0, ExcEmptyObject()); - if (vec_size!=0) - parallel::apply_to_subranges (0U, size(), + 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, _1, _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; } @@ -1126,7 +1129,8 @@ Vector > & Vector >::operator = (const std::complex s) { Assert (numbers::is_finite(s), - ExcMessage("The given value is not finite but either infinite or Not A Number (NaN)")); + ExcMessage("The given value is not finite but either infinite or " + "Not A Number (NaN)")); if (s != std::complex()) Assert (vec_size!=0, ExcEmptyObject()); @@ -1146,13 +1150,15 @@ Vector::operator = (const Vector& v) { if (v.vec_size != vec_size) reinit (v.vec_size, true); - if (vec_size!=0) - parallel::apply_to_subranges (0U, size(), + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::apply_to_subranges (0U, vec_size, std_cxx1x::bind(&internal::Vector::template copy_subrange, std_cxx1x::cref(v), _1, _2, std_cxx1x::ref(*this)), internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + internal::Vector::copy_subrange(v, 0U, vec_size, *this); return *this; } @@ -1165,16 +1171,18 @@ inline Vector & Vector::operator = (const Vector& v) { - if (v.size() != vec_size) - reinit (v.size(), true); - if (vec_size!=0) - parallel::apply_to_subranges (0U, size(), + if (v.vec_size != vec_size) + reinit (v.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_ext, std_cxx1x::cref(v), _1, _2, std_cxx1x::ref(*this)), internal::Vector::minimum_parallel_grain_size); - + else if (vec_size > 0) + internal::Vector::copy_subrange_ext(v, 0U, vec_size, *this); + return *this; } diff --git a/deal.II/lac/include/lac/vector.templates.h b/deal.II/lac/include/lac/vector.templates.h index 9e817471ee..20bf71324e 100644 --- a/deal.II/lac/include/lac/vector.templates.h +++ b/deal.II/lac/include/lac/vector.templates.h @@ -35,6 +35,8 @@ DEAL_II_NAMESPACE_OPEN +#define BLOCK_LEVEL 6 + namespace internal { template @@ -326,16 +328,70 @@ Number Vector::operator * (const Vector& v) const Assert (vec_size == v.size(), ExcDimensionMismatch(vec_size, v.size())); - Number sum = 0; + 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; // 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 - for (unsigned int i=0; i::conjugate(v.val[i])); + // 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) + { + sum1 = 0.; + for (unsigned int j=0; j::conjugate(*Y++)); + sum2 += sum3; + } + sum1 += sum2; + } + sum += sum1; + } + while (X != X_end2) + { + sum2 = 0.; + for (unsigned int j=0; j::conjugate(*Y++)); + sum2 += sum3; + } + sum += sum2; + } + while (X != X_end3) + { + sum3 = 0; + for (unsigned int i=0; i::conjugate(*Y++)); + sum += sum3; + } + + sum3 = 0; + while (X != X_end) + sum3 += *X++ * Number(numbers::NumberTraits::conjugate(*Y++)); + sum += sum3; return sum; } @@ -347,10 +403,54 @@ Vector::norm_sqr () const { Assert (vec_size!=0, ExcEmptyObject()); - real_type sum = 0; + 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)); - for (unsigned int i=0; i::abs_square(val[i]); + while (X != X_end1) + { + sum1 = 0.; + for (unsigned int j=0; j::abs_square(*X++); + sum2 += sum3; + } + sum1 += sum2; + } + sum += sum1; + } + while (X != X_end2) + { + sum2 = 0.; + for (unsigned int j=0; j::abs_square(*X++); + sum2 += sum3; + } + 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; return sum; } @@ -361,10 +461,54 @@ Number Vector::mean_value () const { Assert (vec_size!=0, ExcEmptyObject()); - Number sum = 0; + 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)); - for (unsigned int i=0; i::real_type 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; + } - real_type sum = 0; - - for (unsigned int i=0; i::abs(val[i]); + sum3 = 0; + while (X != X_end) + sum3 += numbers::NumberTraits::abs(*X++); + sum += sum3; return sum; } @@ -402,10 +590,59 @@ Vector::lp_norm (const real_type p) const { Assert (vec_size!=0, ExcEmptyObject()); - real_type sum = 0; + if (p == 1.) + return l1_norm(); + else if (p == 2.) + return std::sqrt(norm_sqr()); - for (unsigned int i=0; i::abs(val[i]), p); + 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++),p); + sum2 += sum3; + } + sum1 += sum2; + } + sum += sum1; + } + while (X != X_end2) + { + sum2 = 0.; + for (unsigned int j=0; j::abs(*X++),p); + sum2 += sum3; + } + sum += sum2; + } + 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)); } @@ -443,12 +680,16 @@ Vector& Vector::operator -= (const Vector& v) 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 - boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + v.val, + val, + (boost::lambda::_1 - boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::add (const Number v) { Assert (vec_size!=0, ExcEmptyObject()); - parallel::transform (val, - val+vec_size, - val, - (boost::lambda::_1 + v), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + val, + (boost::lambda::_1 + v), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::add (const Vector& v) 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 + boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + v.val, + val, + (boost::lambda::_1 + boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::add (const Number a, const Vector& v, Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size)); - parallel::transform (val, - val+vec_size, - v.val, - w.val, - val, - (boost::lambda::_1 + a*boost::lambda::_2 + b*boost::lambda::_3), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + v.val, + w.val, + val, + (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::sadd (const Number x, 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 + boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + v.val, + val, + (x*boost::lambda::_1 + boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::sadd (const Number x, const Number a, Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); 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, @@ -548,6 +806,9 @@ void Vector::sadd (const Number x, const Number a, 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::scale (const Vector &s) Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_size)); - parallel::transform (val, - val+vec_size, - s.val, - val, - (boost::lambda::_1*boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (val, + val+vec_size, + s.val, + val, + (boost::lambda::_1*boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::equ (const Number a, Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size)); - parallel::transform (u.val, - u.val+u.vec_size, - val, - (a*boost::lambda::_1), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (u.val, + u.val+u.vec_size, + val, + (a*boost::lambda::_1), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::equ (const Number a, const Vector& u, Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_size)); Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); - parallel::transform (u.val, - u.val+u.vec_size, - v.val, - val, - (a*boost::lambda::_1 + b*boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); - } + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (u.val, + u.val+u.vec_size, + v.val, + val, + (a*boost::lambda::_1 + b*boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i @@ -666,13 +939,17 @@ void Vector::equ (const Number a, const Vector& u, Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); Assert (vec_size == w.vec_size, ExcDimensionMismatch(vec_size, w.vec_size)); - parallel::transform (u.val, - u.val+u.vec_size, - v.val, - w.val, - val, - (a*boost::lambda::_1 + b*boost::lambda::_2 + c*boost::lambda::_3), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (u.val, + u.val+u.vec_size, + v.val, + w.val, + val, + (a*boost::lambda::_1 + b*boost::lambda::_2 + c*boost::lambda::_3), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i::ratio (const Vector &a, // we overwrite them anyway reinit (a.size(), true); - parallel::transform (a.val, - a.val+a.vec_size, - b.val, - val, - (boost::lambda::_1 / boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + if (vec_size>internal::Vector::minimum_parallel_grain_size) + parallel::transform (a.val, + a.val+a.vec_size, + b.val, + val, + (boost::lambda::_1 / boost::lambda::_2), + internal::Vector::minimum_parallel_grain_size); + else if (vec_size > 0) + for (unsigned int i=0; i