From 78bc2b28106bd1d702eee4dc93858f05d1c3e832 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Wed, 13 Apr 2016 15:42:03 +0200 Subject: [PATCH] Improve grainsize of vectors, unify vector loops --- include/deal.II/lac/vector.templates.h | 167 ++++++++++++------------- source/base/parallel.cc | 39 +++--- 2 files changed, 100 insertions(+), 106 deletions(-) diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 8a3bc24aa5..3bd2faca6e 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -122,10 +122,11 @@ namespace internal void vectorized_transform(Functor &functor, size_type vec_size) { -#ifndef DEAL_II_WITH_THREADS - functor(0,vec_size); -#else - if (vec_size>internal::Vector::minimum_parallel_grain_size) +#ifdef DEAL_II_WITH_THREADS + // only go to the parallel function in case there are at least 4 parallel + // items, otherwise the overhead is too large + if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size && + MultithreadInfo::n_threads() > 1) { tbb::parallel_for (tbb::blocked_range (0, vec_size, @@ -134,12 +135,62 @@ namespace internal tbb::auto_partitioner()); } else if (vec_size > 0) - functor(0,vec_size); #endif + functor(0,vec_size); } - // Define the functors necessary to use SIMD with TBB. + // Define the functors necessary to use SIMD with TBB. we also include the + // simple copy and set operations + + template + struct Vector_set + { + Number *dst; + Number value; + +#ifdef DEAL_II_WITH_THREADS + void operator() (const tbb::blocked_range &range) const + { + operator()(range.begin(),range.end()); + } +#endif + + void operator() (const size_type begin, const size_type end) const + { + if (value == Number()) + std::memset (dst+begin,0,(end-begin)*sizeof(Number)); + else + std::fill (dst+begin, dst+end, value); + } + }; + + template + struct Vector_copy + { + const OtherNumber *src; + Number *dst; + +#ifdef DEAL_II_WITH_THREADS + void operator() (const tbb::blocked_range &range) const + { + operator()(range.begin(),range.end()); + } +#endif + + void operator() (const size_type begin, const size_type end) const + { + if (types_are_equal::value) + std::memcpy(dst+begin, src+begin, (end-begin)*sizeof(Number)); + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (typename dealii::Vector::size_type i=begin; i struct Vectorization_multiply_factor { @@ -554,6 +605,26 @@ namespace internal } }; + // function used in the header file + namespace Vector + { + template + void copy_vector (const dealii::Vector &src, + dealii::Vector &dst) + { + if (PointerComparison::equal(&src, &dst)) + return; + + const typename dealii::Vector::size_type vec_size = src.size(); + const typename dealii::Vector::size_type dst_size = dst.size(); + if (dst_size != vec_size) + dst.reinit (src, true); + dealii::internal::Vector_copy copier; + copier.dst = dst.begin(); + copier.src = src.begin(); + internal::vectorized_transform(copier,vec_size); + } + } } @@ -761,76 +832,6 @@ Vector::is_non_negative () const -namespace internal -{ - namespace Vector - { - template - void set_subrange (const T s, - const typename dealii::Vector::size_type begin, - const typename dealii::Vector::size_type 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 typename dealii::Vector::size_type begin, - const typename dealii::Vector::size_type end, - const dealii::Vector &src, - dealii::Vector &dst) - { - memcpy(&*(dst.begin()+begin), &*(src.begin()+begin), - (end-begin)*sizeof(T)); - } - - - template - void copy_subrange (const typename dealii::Vector::size_type begin, - const typename dealii::Vector::size_type end, - const dealii::Vector &src, - dealii::Vector &dst) - { - const T *src_ptr = src.begin(); - U *dst_ptr = dst.begin(); - DEAL_II_OPENMP_SIMD_PRAGMA - for (typename dealii::Vector::size_type i=begin; i - void copy_vector (const dealii::Vector &src, - dealii::Vector &dst) - { - if (PointerComparison::equal(&src, &dst)) - return; - - const typename dealii::Vector::size_type vec_size = src.size(); - const typename dealii::Vector::size_type 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_cxx11::bind(&internal::Vector::template - copy_subrange , - std_cxx11::_1, - std_cxx11::_2, - std_cxx11::cref(src), - std_cxx11::ref(dst)), - internal::Vector::minimum_parallel_grain_size); - else if (vec_size > 0) - copy_subrange (0U, vec_size, src, dst); - } - } -} - - - template Vector & Vector::operator= (const Number s) @@ -838,14 +839,12 @@ Vector::operator= (const Number s) AssertIsFinite(s); 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_cxx11::bind(&internal::Vector::template - set_subrange, - s, std_cxx11::_1, std_cxx11::_2, std_cxx11::ref(*this)), - internal::Vector::minimum_parallel_grain_size); - else if (vec_size > 0) - internal::Vector::set_subrange(s, 0U, vec_size, *this); + + internal::Vector_set setter; + setter.dst = val; + setter.value = s; + + internal::vectorized_transform(setter,vec_size); return *this; } diff --git a/source/base/parallel.cc b/source/base/parallel.cc index 0ca6d65756..09a175aa9a 100644 --- a/source/base/parallel.cc +++ b/source/base/parallel.cc @@ -23,33 +23,28 @@ namespace internal { namespace Vector { - // set minimum grain size. this value is - // roughly in accordance with the curve - // in the TBB book (fig 3.2) that shows - // run time as a function of grain size - // -- there, values from 200 upward are - // so that the scheduling overhead - // amortizes well (for very large values - // in that example, the grain size is too - // large to split the work load into - // enough chunks and the problem becomes - // badly balanced) - unsigned int minimum_parallel_grain_size = 1024; + // set minimum grain size. this value has been determined by experiments + // with the actual dealii::Vector implementation and takes vectorization + // that is done inside most functions of dealii::Vector into account + // (without vectorization, we would land at around 1000). for smaller + // values than this, the scheduling overhead will be significant. if the + // value becomes too large, on the other hand, the work load cannot be + // split into enough chunks and the problem becomes badly balanced. the + // tests are documented at https://github.com/dealii/dealii/issues/2496 + unsigned int minimum_parallel_grain_size = 4096; } namespace SparseMatrix { - // set this value to 1/4 of the value of - // the minimum grain size of - // vectors. this rests on the fact that - // we have to do a lot more work per row - // of a matrix than per element of a - // vector. it could possibly be reduced - // even further but that doesn't appear - // worth it any more for anything but - // very small matrices that we don't care - // that much about anyway. + // set this value to 1/16 of the value of the minimum grain size of + // vectors (a factor of 4 because we assume that sparse matrix-vector + // products do not vectorize and the other factor of 4 because we expect + // at least four entries per row). this rests on the fact that we have to + // do a lot more work per row of a matrix than per element of a vector. it + // could possibly be reduced even further but that doesn't appear worth it + // any more for anything but very small matrices that we don't care that + // much about anyway. unsigned int minimum_parallel_grain_size = 256; } } -- 2.39.5