From 9cf406844e67c2d545acc4c9588ff2289a1ab6a6 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Thu, 14 Apr 2016 16:16:13 +0200 Subject: [PATCH] Implement affinity partitioner for the deal.II vector class --- include/deal.II/base/parallel.h | 73 ++++++++++ include/deal.II/lac/vector.h | 99 ++----------- include/deal.II/lac/vector.templates.h | 188 +++++++++++++++++++------ source/lac/vector.cc | 17 +-- source/lac/vector.inst.in | 27 +--- 5 files changed, 240 insertions(+), 164 deletions(-) diff --git a/include/deal.II/base/parallel.h b/include/deal.II/base/parallel.h index 164afa8272..2aeab248b3 100644 --- a/include/deal.II/base/parallel.h +++ b/include/deal.II/base/parallel.h @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -683,6 +684,78 @@ namespace parallel #endif } + +// --------------------- for loop affinity partitioner ----------------------- + + /** + * A class that wraps a TBB affinity partitioner in a thread-safe way. In + * Vector, we use a shared pointer to share an affinity partitioner + * between different vectors of the same size for improving data (and + * NUMA) locality. However, when an outer task does multiple vector + * operations, the shared pointer could lead to race conditions. This + * class only allows one instance to get a partitioner. The other objects + * cannot use that object and need to create their own copy. + */ + namespace internal + { + class TBBPartitioner + { + public: + /** + * Constructor. + */ + TBBPartitioner() +#ifdef DEAL_II_WITH_THREADS + : + in_use(false) +#endif + {} + +#ifdef DEAL_II_WITH_THREADS + /** + * Return a pointer to the affinity partitioner in case it is free. If + * it is already in use (i.e., someone else called acquire without + * releasing it afterwards), the null pointer is returned. + */ + tbb::affinity_partitioner *acquire() + { + dealii::Threads::Mutex::ScopedLock lock(mutex); + if (in_use) + return NULL; + + in_use = true; + return &partitioner; + } + + /** + * After using the partitioner in a tbb loop through acquire(), this + * call makes the partitioner available again. + */ + void release() + { + dealii::Threads::Mutex::ScopedLock lock(mutex); + in_use = false; + } + + private: + /** + * The stored partitioner + */ + tbb::affinity_partitioner partitioner; + + /** + * A flag to indicate whether the partitioner has been acquired but not + * released yet, i.e., it is in use somewhere else. + */ + bool in_use; + + /** + * A mutex to guard the access to the in_use flag. + */ + dealii::Threads::Mutex mutex; +#endif + }; + } } diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index e507c90ca8..fea7f5a1ec 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -61,6 +61,13 @@ template class BlockVector; template class VectorView; +namespace parallel +{ + namespace internal + { + class TBBPartitioner; + } +} @@ -984,6 +991,12 @@ protected: */ Number *val; + /** + * For parallel loops with TBB, this member variable stores the affinity + * information of loops. + */ + std_cxx11::shared_ptr thread_loop_partitioner; + /** * Make all other vector types friends. */ @@ -1073,92 +1086,6 @@ Vector::~Vector () -template -inline -void Vector::reinit (const size_type n, - const bool omit_zeroing_entries) -{ - if (n==0) - { - if (val) deallocate(); - val = 0; - max_vec_size = vec_size = 0; - return; - }; - - if (n>max_vec_size) - { - if (val) deallocate(); - max_vec_size = n; - allocate(); - }; - vec_size = n; - if (omit_zeroing_entries == false) - *this = static_cast(0); -} - - - -// declare function that is implemented in vector.templates.h -namespace internal -{ - namespace Vector - { - template - void copy_vector (const dealii::Vector &src, - dealii::Vector &dst); - } -} - - - -template -inline -Vector & -Vector::operator= (const Vector &v) -{ - dealii::internal::Vector::copy_vector (v, *this); - return *this; -} - - - -#ifdef DEAL_II_WITH_CXX11 -template -inline -Vector & -Vector::operator= (Vector &&v) -{ - Subscriptor::operator=(std::move(v)); - - if (val) deallocate(); - - vec_size = v.vec_size; - max_vec_size = v.max_vec_size; - val = v.val; - - v.vec_size = 0; - v.max_vec_size = 0; - v.val = nullptr; - - return *this; -} -#endif - - - -template -template -inline -Vector & -Vector::operator= (const Vector &v) -{ - internal::Vector::copy_vector (v, *this); - return *this; -} - - - template inline std::size_t Vector::size () const diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index 24f0416fa9..aff70219c6 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -120,7 +120,8 @@ namespace internal template void vectorized_transform(Functor &functor, - size_type vec_size) + size_type vec_size, + parallel::internal::TBBPartitioner *partitioner = NULL) { #ifdef DEAL_II_WITH_THREADS // only go to the parallel function in case there are at least 4 parallel @@ -128,15 +129,29 @@ namespace internal if (vec_size >= 4*internal::Vector::minimum_parallel_grain_size && MultithreadInfo::n_threads() > 1) { + // in case the partitioner was used by another thread, create a + // private partitioner here for use in the loop (we could also use the + // auto partitioner of TBB but that increases code size considerably + // without major benefit) + tbb::affinity_partitioner *tbb_partitioner = (partitioner == NULL) ? + NULL : partitioner->acquire(); + tbb::affinity_partitioner private_partitioner; + if (tbb_partitioner == NULL) + tbb_partitioner = &private_partitioner; + tbb::parallel_for (tbb::blocked_range (0, vec_size, internal::Vector::minimum_parallel_grain_size), functor, - tbb::auto_partitioner()); + *tbb_partitioner); + if (partitioner != NULL) + partitioner->release(); } else if (vec_size > 0) -#endif functor(0,vec_size); +#else + functor(0,vec_size); +#endif } @@ -604,27 +619,6 @@ 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); - } - } } @@ -653,7 +647,8 @@ Vector::Vector (Vector &&v) Subscriptor(std::move(v)), vec_size(v.vec_size), max_vec_size(v.max_vec_size), - val(v.val) + val(v.val), + thread_loop_partitioner(std::move(v.thread_loop_partitioner)) { v.vec_size = 0; v.max_vec_size = 0; @@ -795,12 +790,123 @@ Vector::Vector (const TrilinosWrappers::Vector &v) #endif +template +inline +Vector & +Vector::operator= (const Vector &v) +{ + if (PointerComparison::equal(this, &v)) + return *this; + + if (vec_size != v.vec_size) + reinit (v, true); + + dealii::internal::Vector_copy copier; + copier.dst = val; + copier.src = v.val; + internal::vectorized_transform(copier,vec_size,thread_loop_partitioner.get()); + + return *this; +} + + + +#ifdef DEAL_II_WITH_CXX11 +template +inline +Vector & +Vector::operator= (Vector &&v) +{ + Subscriptor::operator=(std::move(v)); + + if (val) deallocate(); + + vec_size = v.vec_size; + max_vec_size = v.max_vec_size; + val = v.val; + thread_loop_partitioner = std::move(v.thread_loop_partitioner); + + v.vec_size = 0; + v.max_vec_size = 0; + v.val = nullptr; + + return *this; +} +#endif + + + +template +template +inline +Vector & +Vector::operator= (const Vector &v) +{ + if (vec_size != v.vec_size) + reinit (v, true); + + dealii::internal::Vector_copy copier; + copier.dst = val; + copier.src = v.val; + internal::vectorized_transform(copier,vec_size,thread_loop_partitioner.get()); + + return *this; +} + + + +template +inline +void Vector::reinit (const size_type n, + const bool omit_zeroing_entries) +{ + thread_loop_partitioner.reset(); + if (n==0) + { + if (val) deallocate(); + val = 0; + max_vec_size = vec_size = 0; + return; + }; + + if (n>max_vec_size) + { + if (val) deallocate(); + max_vec_size = n; + allocate(); + }; + vec_size = n; + thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner()); + if (omit_zeroing_entries == false) + *this = static_cast(0); +} + + + template template void Vector::reinit (const Vector &v, const bool omit_zeroing_entries) { - reinit (v.size(), omit_zeroing_entries); + thread_loop_partitioner = v.thread_loop_partitioner; + + if (v.vec_size==0) + { + if (val) deallocate(); + val = 0; + max_vec_size = vec_size = 0; + return; + }; + + if (v.vec_size>max_vec_size) + { + if (val) deallocate(); + max_vec_size = v.vec_size; + allocate(); + }; + vec_size = v.vec_size; + if (omit_zeroing_entries == false) + *this = static_cast(0); } @@ -846,7 +952,7 @@ Vector::operator= (const Number s) setter.dst = val; setter.value = s; - internal::vectorized_transform(setter,vec_size); + internal::vectorized_transform(setter,vec_size,thread_loop_partitioner.get()); return *this; } @@ -881,7 +987,7 @@ Vector &Vector::operator *= (const Number factor) vector_multiply.val = val; vector_multiply.factor = factor; - internal::vectorized_transform(vector_multiply,vec_size); + internal::vectorized_transform(vector_multiply,vec_size,thread_loop_partitioner.get()); return *this; } @@ -902,7 +1008,7 @@ Vector::add (const Number a, vector_add_av.val = val; vector_add_av.v_val = v.val; vector_add_av.factor = a; - internal::vectorized_transform(vector_add_av,vec_size); + internal::vectorized_transform(vector_add_av,vec_size,thread_loop_partitioner.get()); } @@ -924,7 +1030,7 @@ Vector::sadd (const Number x, vector_sadd_xav.v_val = v.val; vector_sadd_xav.a = a; vector_sadd_xav.x = x; - internal::vectorized_transform(vector_sadd_xav,vec_size); + internal::vectorized_transform(vector_sadd_xav,vec_size,thread_loop_partitioner.get()); } @@ -1540,7 +1646,7 @@ Vector &Vector::operator -= (const Vector &v) internal::Vectorization_subtract_v vector_subtract; vector_subtract.val = val; vector_subtract.v_val = v.val; - internal::vectorized_transform(vector_subtract,vec_size); + internal::vectorized_transform(vector_subtract,vec_size,thread_loop_partitioner.get()); return *this; } @@ -1555,7 +1661,7 @@ void Vector::add (const Number v) internal::Vectorization_add_factor vector_add; vector_add.val = val; vector_add.factor = v; - internal::vectorized_transform(vector_add,vec_size); + internal::vectorized_transform(vector_add,vec_size,thread_loop_partitioner.get()); } @@ -1569,7 +1675,7 @@ void Vector::add (const Vector &v) internal::Vectorization_add_v vector_add; vector_add.val = val; vector_add.v_val = v.val; - internal::vectorized_transform(vector_add,vec_size); + internal::vectorized_transform(vector_add,vec_size,thread_loop_partitioner.get()); } @@ -1591,7 +1697,7 @@ void Vector::add (const Number a, const Vector &v, vector_add.w_val = w.val; vector_add.a = a; vector_add.b = b; - internal::vectorized_transform(vector_add,vec_size); + internal::vectorized_transform(vector_add,vec_size,thread_loop_partitioner.get()); } @@ -1609,7 +1715,7 @@ void Vector::sadd (const Number x, vector_sadd.val = val; vector_sadd.v_val = v.val; vector_sadd.x = x; - internal::vectorized_transform(vector_sadd,vec_size); + internal::vectorized_transform(vector_sadd,vec_size,thread_loop_partitioner.get()); } @@ -1634,7 +1740,7 @@ void Vector::sadd (const Number x, const Number a, vector_sadd.x = x; vector_sadd.a = a; vector_sadd.b = b; - internal::vectorized_transform(vector_sadd,vec_size); + internal::vectorized_transform(vector_sadd,vec_size,thread_loop_partitioner.get()); } @@ -1659,7 +1765,7 @@ void Vector::scale (const Vector &s) internal::Vectorization_scale vector_scale; vector_scale.val = val; vector_scale.v_val = s.val; - internal::vectorized_transform(vector_scale,vec_size); + internal::vectorized_transform(vector_scale,vec_size,thread_loop_partitioner.get()); } @@ -1690,7 +1796,7 @@ void Vector::equ (const Number a, vector_equ.val = val; vector_equ.u_val = u.val; vector_equ.a = a; - internal::vectorized_transform(vector_equ,vec_size); + internal::vectorized_transform(vector_equ,vec_size,thread_loop_partitioner.get()); } @@ -1734,7 +1840,7 @@ void Vector::equ (const Number a, const Vector &u, vector_equ.v_val = v.val; vector_equ.a = a; vector_equ.b = b; - internal::vectorized_transform(vector_equ,vec_size); + internal::vectorized_transform(vector_equ,vec_size,thread_loop_partitioner.get()); } @@ -1756,7 +1862,7 @@ void Vector::equ (const Number a, const Vector &u, vector_equ.a = a; vector_equ.b = b; vector_equ.c = c; - internal::vectorized_transform(vector_equ,vec_size); + internal::vectorized_transform(vector_equ,vec_size,thread_loop_partitioner.get()); } @@ -1776,7 +1882,7 @@ void Vector::ratio (const Vector &a, vector_ratio.val = val; vector_ratio.a_val = a.val; vector_ratio.b_val = b.val; - internal::vectorized_transform(vector_ratio,vec_size); + internal::vectorized_transform(vector_ratio,vec_size,thread_loop_partitioner.get()); } diff --git a/source/lac/vector.cc b/source/lac/vector.cc index 9da88f6780..e7bf7fbe34 100644 --- a/source/lac/vector.cc +++ b/source/lac/vector.cc @@ -21,17 +21,7 @@ DEAL_II_NAMESPACE_OPEN // instantiate for integers: template class Vector; -namespace internal -{ - namespace Vector - { - template void copy_vector (const dealii::Vector &, - dealii::Vector &); - - template void copy_vector (const dealii::Vector &, - dealii::Vector &); - } -} +template Vector &Vector::operator= (const dealii::Vector &); template void Vector::reinit(const Vector &, const bool); @@ -42,8 +32,9 @@ void Vector::reinit(const Vector &, const bool); // arguments is covered by the default copy constructor and copy operator that // is declared separately) -#define TEMPL_COPY_CONSTRUCTOR(S1,S2) \ - template Vector::Vector (const Vector &) +#define TEMPL_COPY_CONSTRUCTOR(S1,S2) \ + template Vector::Vector (const Vector &); \ + template Vector& Vector::operator= (const Vector &) #ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG TEMPL_COPY_CONSTRUCTOR(double,float); diff --git a/source/lac/vector.inst.in b/source/lac/vector.inst.in index 4a7b11874b..34078f5d90 100644 --- a/source/lac/vector.inst.in +++ b/source/lac/vector.inst.in @@ -22,14 +22,6 @@ 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; @@ -49,14 +41,6 @@ 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; @@ -69,11 +53,6 @@ for (S1, S2 : COMPLEX_SCALARS) for (S1: REAL_SCALARS; S2: COMPLEX_SCALARS) { - namespace internal - \{ - namespace Vector - \{ - template void copy_vector (const dealii::Vector&, - dealii::Vector&); - \} - \} + template + Vector& Vector::operator= (const Vector &); + } \ No newline at end of file -- 2.39.5