From 93bd8983076ce65f9c7a78e25645e983b1fc09b5 Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Thu, 21 Aug 2014 18:48:57 -0500 Subject: [PATCH] Replace transform by vectorized_transform. Replace initialize_val by allocate and clear_val by deallocate. Add vectorization to changes.h --- cmake/setup_compiler_flags_gnu.cmake | 2 +- cmake/setup_compiler_flags_intel.cmake | 2 +- doc/news/changes.h | 6 + include/deal.II/lac/vector.h | 50 +-- include/deal.II/lac/vector.templates.h | 596 +++++++++++++++++-------- 5 files changed, 434 insertions(+), 222 deletions(-) diff --git a/cmake/setup_compiler_flags_gnu.cmake b/cmake/setup_compiler_flags_gnu.cmake index f5017fd2cf..21159946c1 100644 --- a/cmake/setup_compiler_flags_gnu.cmake +++ b/cmake/setup_compiler_flags_gnu.cmake @@ -91,7 +91,7 @@ ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-deprecated-declarations") # warning, if OpenMP 4.0 is not supported. # ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-fopenmp-simd") -ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unkown-pragmas") +ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unknown-pragmas") IF(CMAKE_CXX_COMPILER_ID MATCHES "Clang") # diff --git a/cmake/setup_compiler_flags_intel.cmake b/cmake/setup_compiler_flags_intel.cmake index 4b8080cd5e..6961ee6eaa 100644 --- a/cmake/setup_compiler_flags_intel.cmake +++ b/cmake/setup_compiler_flags_intel.cmake @@ -144,7 +144,7 @@ ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-wd280") # warning, if OpenMP 4.0 is not supported. # ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-openmp-simd") -ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unkown-pragmas") +ENABLE_IF_SUPPORTED(DEAL_II_CXX_FLAGS "-Wno-unknown-pragmas") diff --git a/doc/news/changes.h b/doc/news/changes.h index a2151b56d0..51faa4b861 100644 --- a/doc/news/changes.h +++ b/doc/news/changes.h @@ -97,6 +97,12 @@ inconvenience this causes.
    +
  1. New: Most of the operations done on Vector (like add, sadd, etc.) + are now vectorized (SIMD) using OpenMP 4.0 +
    + (Bruno Turcksin, 2014/08/25) +
    2. +
  2. Updated: The version of BOOST that comes bundled with deal.II has been updated to 1.56.0.
    diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index 36076effb1..5e6cf301a3 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -31,10 +31,6 @@ #include #include -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 -#include -#endif - DEAL_II_NAMESPACE_OPEN @@ -962,14 +958,14 @@ protected: private: /** - * Allocate @p v and, if possible, align along 64-byte boundaries. + * Allocate and align @p v along 64-byte boundaries. */ - void initialize_val(const size_type n); + void allocate(const size_type n); /** - * Free @p val. + * Deallocate @p val. */ - void clear_val(); + void deallocate(); }; /*@}*/ @@ -1025,7 +1021,7 @@ Vector::~Vector () { if (val) { - clear_val(); + deallocate(); val=0; } } @@ -1038,7 +1034,7 @@ void Vector::reinit (const size_type n, const bool fast) { if (n==0) { - if (val) clear_val(); + if (val) deallocate(); val = 0; max_vec_size = vec_size = 0; return; @@ -1046,8 +1042,8 @@ void Vector::reinit (const size_type n, const bool fast) if (n>max_vec_size) { - if (val) clear_val(); - initialize_val(n); + if (val) deallocate(); + allocate(n); Assert (val != 0, ExcOutOfMemory()); max_vec_size = n; }; @@ -1349,38 +1345,10 @@ Vector::load (Archive &ar, const unsigned int) ar &vec_size &max_vec_size ; - initialize_val(max_vec_size); + allocate(max_vec_size); ar &boost::serialization::make_array(val, max_vec_size); } - - -template -inline -void -Vector::initialize_val(const size_type size) -{ -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - val = static_cast(_mm_malloc (size, 64)); -#else - val = static_cast(malloc (size)); -#endif -} - - - -template -inline -void -Vector::clear_val() -{ -#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 - _mm_free(val); -#else - free(val); -#endif -} - #endif diff --git a/include/deal.II/lac/vector.templates.h b/include/deal.II/lac/vector.templates.h index ec034a507e..80275ae3f2 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -35,12 +35,13 @@ # include #endif -#include + #include #include #include #include #include +#include DEAL_II_NAMESPACE_OPEN @@ -118,21 +119,330 @@ namespace internal "into a vector of reals/doubles")); } + template + 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) + { + tbb::parallel_for (tbb::blocked_range (0, + vec_size, + internal::Vector::minimum_parallel_grain_size), + functor, + tbb::auto_partitioner()); + } + else if (vec_size > 0) + functor(0,vec_size); +#endif + } + + + // Define the functors neccessary to use SIMD with TBB. + template + struct Vectorization_multiply_factor + { + Number *val; + Number factor; + +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_add_av + { + Number *val; + Number *v_val; + Number factor; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_sadd_xav + { + Number *val; + Number *v_val; + Number a; + Number x; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_subtract_v + { + Number *val; + Number *v_val; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_add_factor + { + Number *val; + Number factor; #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 + { + #pragma omp simd + for (size_type i=begin; i - struct Vectorization_add_1 + struct Vectorization_add_v { - Number* val; - Number* v_val; + Number *val; + Number *v_val; +#ifdef DEAL_II_WITH_THREADS void operator() (const tbb::blocked_range &range) const { -#pragma omp simd - for (size_type i=range.begin(); i!=range.end(); ++i) + operator()(range.begin(),range.end()); + } +#endif + + void operator() (const size_type begin, const size_type end) const + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_add_avpbw + { + Number *val; + Number *v_val; + Number *w_val; + Number a; + Number b; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_sadd_xv + { + Number *val; + Number *v_val; + Number x; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_sadd_xavbw + { + Number *val; + Number *v_val; + Number *w_val; + Number x; + Number a; + Number b; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_scale + { + Number *val; + Number *v_val; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_equ_au + { + Number *val; + Number *u_val; + Number a; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_equ_aubv + { + Number *val; + Number *u_val; + Number *v_val; + Number a; + Number b; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_equ_aubvcw + { + Number *val; + Number *u_val; + Number *v_val; + Number *w_val; + Number a; + Number b; + Number c; +#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 + { + #pragma omp simd + for (size_type i=begin; i + struct Vectorization_ratio + { + Number *val; + Number *a_val; + Number *b_val; +#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 + { + #pragma omp simd + for (size_type i=begin; i::Vector (const Vector &v) { if (vec_size != 0) { - initialize_val(max_vec_size); + allocate(max_vec_size); Assert (val != 0, ExcOutOfMemory()); *this = v; } @@ -168,7 +478,7 @@ Vector::Vector (const Vector &v) { if (vec_size != 0) { - initialize_val(max_vec_size); + allocate(max_vec_size); Assert (val != 0, ExcOutOfMemory()); std::copy (v.begin(), v.end(), begin()); } @@ -241,7 +551,7 @@ Vector::Vector (const TrilinosWrappers::MPI::Vector &v) { if (vec_size != 0) { - initialize_val(max_vec_size); + allocate(max_vec_size); Assert (val != 0, ExcOutOfMemory()); // Copy the distributed vector to @@ -275,7 +585,7 @@ Vector::Vector (const TrilinosWrappers::Vector &v) { if (vec_size != 0) { - initialize_val(max_vec_size); + allocate(max_vec_size); Assert (val != 0, ExcOutOfMemory()); // get a representation of the vector @@ -384,16 +694,6 @@ namespace internal } - template - void copy_subrange_wrap (const typename dealii::Vector::size_type begin, - const typename dealii::Vector::size_type end, - const dealii::Vector &src, - dealii::Vector &dst) - { - copy_subrange (begin, end, src, dst); - } - - template void copy_vector (const dealii::Vector &src, dealii::Vector &dst) @@ -408,7 +708,7 @@ namespace internal 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 , + copy_subrange , std_cxx1x::_1, std_cxx1x::_2, std_cxx1x::cref(src), @@ -467,11 +767,11 @@ Vector &Vector::operator *= (const Number factor) Assert (vec_size!=0, ExcEmptyObject()); - parallel::transform (val, - val+vec_size, - val, - (factor*boost::lambda::_1), - internal::Vector::minimum_parallel_grain_size); + internal::Vectorization_multiply_factor vector_multiply; + vector_multiply.val = val; + vector_multiply.factor = factor; + + internal::vectorized_transform(vector_multiply,vec_size); return *this; } @@ -488,12 +788,11 @@ Vector::add (const Number a, 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); + internal::Vectorization_add_av vector_add_av; + 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); } @@ -510,12 +809,12 @@ Vector::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 + a*boost::lambda::_2), - internal::Vector::minimum_parallel_grain_size); + internal::Vectorization_sadd_xav vector_sadd_xav; + vector_sadd_xav.val = val; + 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); } @@ -934,16 +1233,10 @@ Vector &Vector::operator -= (const Vector &v) Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_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 (size_type i=0; i vector_subtract; + vector_subtract.val = val; + vector_subtract.v_val = v.val; + internal::vectorized_transform(vector_subtract,vec_size); return *this; } @@ -954,15 +1247,10 @@ void Vector::add (const Number v) { Assert (vec_size!=0, ExcEmptyObject()); - 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 (size_type i=0; i vector_add; + vector_add.val = val; + vector_add.factor = v; + internal::vectorized_transform(vector_add,vec_size); } @@ -972,45 +1260,10 @@ void Vector::add (const Vector &v) Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); - // Vectorization can only be if v is not the current vector otherwise there - // will be collision - if (&(v.val[0])!=&val[0]) - { -#ifndef DEAL_II_WITH_THREADS -#pragma omp simd - for (size_type i=0; iinternal::Vector::minimum_parallel_grain_size) - { - internal::Vectorization_add_1 vector_add; - vector_add.val = val; - vector_add.v_val = v.val; - tbb::parallel_for (tbb::blocked_range (0, - vec_size, - internal::Vector::minimum_parallel_grain_size), - vector_add, - tbb::auto_partitioner()); - } - else if (vec_size > 0) -#pragma omp simd - for (size_type i=0; iinternal::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 (size_type i=0; i vector_add; + vector_add.val = val; + vector_add.v_val = v.val; + internal::vectorized_transform(vector_add,vec_size); } @@ -1025,17 +1278,13 @@ void Vector::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)); - 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 (size_type i=0; i vector_add; + vector_add.val = val; + vector_add.v_val = v.val; + vector_add.w_val = w.val; + vector_add.a = a; + vector_add.b = b; + internal::vectorized_transform(vector_add,vec_size); } @@ -1049,16 +1298,11 @@ void Vector::sadd (const Number x, Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_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 (size_type i=0; i vector_sadd; + vector_sadd.val = val; + vector_sadd.v_val = v.val; + vector_sadd.x = x; + internal::vectorized_transform(vector_sadd,vec_size); } @@ -1076,17 +1320,14 @@ void Vector::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, - 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 (size_type i=0; i vector_sadd; + vector_sadd.val = val; + vector_sadd.v_val = v.val; + vector_sadd.w_val = w.val; + vector_sadd.x = x; + vector_sadd.a = a; + vector_sadd.b = b; + internal::vectorized_transform(vector_sadd,vec_size); } @@ -1108,16 +1349,10 @@ void Vector::scale (const Vector &s) Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == s.vec_size, ExcDimensionMismatch(vec_size, s.vec_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 (size_type i=0; i vector_scale; + vector_scale.val = val; + vector_scale.v_val = s.val; + internal::vectorized_transform(vector_scale,vec_size); } @@ -1144,15 +1379,11 @@ void Vector::equ (const Number a, Assert (vec_size!=0, ExcEmptyObject()); Assert (vec_size == u.vec_size, ExcDimensionMismatch(vec_size, u.vec_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 (size_type i=0; i vector_equ; + vector_equ.val = val; + vector_equ.u_val = u.val; + vector_equ.a = a; + internal::vectorized_transform(vector_equ,vec_size); } @@ -1190,16 +1421,13 @@ void Vector::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)); - 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 (size_type i=0; i vector_equ; + vector_equ.val = val; + vector_equ.u_val = u.val; + vector_equ.v_val = v.val; + vector_equ.a = a; + vector_equ.b = b; + internal::vectorized_transform(vector_equ,vec_size); } @@ -1213,17 +1441,15 @@ 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)); - 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 (size_type i=0; i vector_equ; + vector_equ.val = val; + vector_equ.u_val = u.val; + vector_equ.v_val = v.val; + vector_equ.w_val = w.val; + vector_equ.a = a; + vector_equ.b = b; + vector_equ.c = c; + internal::vectorized_transform(vector_equ,vec_size); } @@ -1239,16 +1465,11 @@ void Vector::ratio (const Vector &a, // we overwrite them anyway reinit (a.size(), true); - 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 (size_type i=0; i vector_ratio; + vector_ratio.val = val; + vector_ratio.a_val = a.val; + vector_ratio.b_val = b.val; + internal::vectorized_transform(vector_ratio,vec_size); } @@ -1524,6 +1745,23 @@ Vector::memory_consumption () const } + +template +void +Vector::allocate(const size_type size) +{ + val = static_cast(_mm_malloc (sizeof(Number)*size, 64)); +} + + + +template +void +Vector::deallocate() +{ + _mm_free(val); +} + DEAL_II_NAMESPACE_CLOSE #endif -- 2.39.5