From d97be5f007e72126c40165376c9dd23419aeff4c Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Thu, 21 Aug 2014 10:12:54 -0500 Subject: [PATCH] Align the underlying array of Vector and add runtime switched between add and vectorized add. --- include/deal.II/lac/vector.h | 61 +++++++++++++++---- include/deal.II/lac/vector.templates.h | 84 +++++++++++++------------- 2 files changed, 90 insertions(+), 55 deletions(-) diff --git a/include/deal.II/lac/vector.h b/include/deal.II/lac/vector.h index 59205f0dc2..36076effb1 100644 --- a/include/deal.II/lac/vector.h +++ b/include/deal.II/lac/vector.h @@ -31,6 +31,10 @@ #include #include +#if DEAL_II_COMPILER_VECTORIZATION_LEVEL > 0 +#include +#endif + DEAL_II_NAMESPACE_OPEN @@ -626,14 +630,6 @@ public: */ void add (const Vector &V); - /** - * Simple vectorized (SIMD) vector addition, equal to the operator +=. - * This function should not be used to add the vector to itself, i.e to - * multiply the current vector by two. - * - * @dealiiOperationIsMultithreaded - */ - void vectorized_add (const Vector &V); /** * Multiple addition of scaled vectors, i.e. *this += a*V+b*W. @@ -962,6 +958,18 @@ protected: * VectorView will access the pointer. */ friend class VectorView; + +private: + + /** + * Allocate @p v and, if possible, align along 64-byte boundaries. + */ + void initialize_val(const size_type n); + + /** + * Free @p val. + */ + void clear_val(); }; /*@}*/ @@ -1017,7 +1025,7 @@ Vector::~Vector () { if (val) { - delete[] val; + clear_val(); val=0; } } @@ -1030,7 +1038,7 @@ void Vector::reinit (const size_type n, const bool fast) { if (n==0) { - if (val) delete[] val; + if (val) clear_val(); val = 0; max_vec_size = vec_size = 0; return; @@ -1038,8 +1046,8 @@ void Vector::reinit (const size_type n, const bool fast) if (n>max_vec_size) { - if (val) delete[] val; - val = new value_type[n]; + if (val) clear_val(); + initialize_val(n); Assert (val != 0, ExcOutOfMemory()); max_vec_size = n; }; @@ -1341,11 +1349,38 @@ Vector::load (Archive &ar, const unsigned int) ar &vec_size &max_vec_size ; - val = new Number[max_vec_size]; + initialize_val(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 0289f655b7..ec034a507e 100644 --- a/include/deal.II/lac/vector.templates.h +++ b/include/deal.II/lac/vector.templates.h @@ -148,7 +148,7 @@ Vector::Vector (const Vector &v) { if (vec_size != 0) { - val = new Number[max_vec_size]; + initialize_val(max_vec_size); Assert (val != 0, ExcOutOfMemory()); *this = v; } @@ -168,7 +168,7 @@ Vector::Vector (const Vector &v) { if (vec_size != 0) { - val = new Number[max_vec_size]; + initialize_val(max_vec_size); Assert (val != 0, ExcOutOfMemory()); std::copy (v.begin(), v.end(), begin()); } @@ -189,7 +189,7 @@ Vector::Vector (const PETScWrappers::Vector &v) { if (vec_size != 0) { - val = new Number[max_vec_size]; + iniatilize_val(max_vec_size); Assert (val != 0, ExcOutOfMemory()); // get a representation of the vector @@ -241,7 +241,7 @@ Vector::Vector (const TrilinosWrappers::MPI::Vector &v) { if (vec_size != 0) { - val = new Number[max_vec_size]; + initialize_val(max_vec_size); Assert (val != 0, ExcOutOfMemory()); // Copy the distributed vector to @@ -275,7 +275,7 @@ Vector::Vector (const TrilinosWrappers::Vector &v) { if (vec_size != 0) { - val = new Number[max_vec_size]; + initialize_val(max_vec_size); Assert (val != 0, ExcOutOfMemory()); // get a representation of the vector @@ -972,16 +972,45 @@ void Vector::add (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) + // 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::add (const Number a, const Vector &v, } -template -void Vector::vectorized_add (const Vector &v) -{ - Assert (vec_size!=0, ExcEmptyObject()); - Assert (vec_size == v.vec_size, ExcDimensionMismatch(vec_size, v.vec_size)); - -#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; i void Vector::sadd (const Number x, -- 2.39.5