From 8a5fe14def02eb9b4e1ea805a17f552549fa2cba Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Wed, 4 May 2016 09:40:27 -0400 Subject: [PATCH] Use multithreading for ReadWriteVector and LinearAlgebra::Vector --- include/deal.II/lac/la_vector.h | 73 +---- include/deal.II/lac/la_vector.templates.h | 258 +++++++++++------- include/deal.II/lac/read_write_vector.h | 71 ++--- .../deal.II/lac/read_write_vector.templates.h | 77 +++++- include/deal.II/lac/vector.templates.h | 13 - include/deal.II/lac/vector_internal.h | 7 +- source/lac/la_vector.cc | 20 ++ source/lac/read_write_vector.cc | 23 ++ tests/lac/readwritevector_assignment.cc | 67 ----- tests/lac/readwritevector_assignment.output | 55 ---- 10 files changed, 307 insertions(+), 357 deletions(-) delete mode 100644 tests/lac/readwritevector_assignment.cc delete mode 100644 tests/lac/readwritevector_assignment.output diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h index c89eb06094..51b1dce502 100644 --- a/include/deal.II/lac/la_vector.h +++ b/include/deal.II/lac/la_vector.h @@ -267,23 +267,6 @@ namespace LinearAlgebra DeclException0(ExcVectorTypeNotCompatible); private: - /** - * Compute the L1 norm in a recursive way by dividing the vector on - * smaller and smaller intervals. This reduces the numerical error on - * large vector. - */ - typename VectorSpaceVector::real_type l1_norm_recursive(unsigned int i, - unsigned int j) const; - - /** - * Compute the squared L2 norm in a recursive way by dividing the vector - * on smaller and smaller intervals. This reduces the numerical error on - * large vector. - */ - typename VectorSpaceVector::real_type l2_norm_squared_recursive( - unsigned int i, - unsigned int j) const; - /** * Serialize the data of this object using boost. This function is * necessary to use boost::archive::text_iarchive and @@ -293,6 +276,11 @@ namespace LinearAlgebra void serialize(Archive &ar, const unsigned int version); friend class boost::serialization::access; + + /** + * Make all other ReadWriteVector types friends. + */ + template friend class Vector; }; /*@}*/ @@ -339,57 +327,6 @@ namespace LinearAlgebra - template - inline - Vector &Vector::operator= (const Vector &in_vector) - { - this->reinit(in_vector.size(),true); - std::copy(in_vector.begin(), in_vector.end(), this->begin()); - - return *this; - } - - - - template - template - inline - Vector &Vector::operator= (const Vector &in_vector) - { - this->reinit(in_vector.size(in_vector.size()),true); - std::copy(in_vector.begin(), in_vector.end(), this->begin()); - - return *this; - } - - - - template - inline - Vector &Vector::operator= (const Number s) - { - Assert(s==static_cast(0), ExcMessage("Only 0 can be assigned to a vector.")); - (void) s; - - std::fill(this->begin(),this->end(),Number()); - - return *this; - } - - - - template - inline - void Vector::add(const Number a) - { - AssertIsFinite(a); - const size_type size = this->size(); - for (size_type i=0; ival[i] += a; - } - - - template inline typename Vector::size_type Vector::size() const diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h index 5743151a64..4ec2db73e8 100644 --- a/include/deal.II/lac/la_vector.templates.h +++ b/include/deal.II/lac/la_vector.templates.h @@ -17,6 +17,7 @@ #define dealii__la_vector_templates_h #include +#include #include #include @@ -24,13 +25,71 @@ DEAL_II_NAMESPACE_OPEN namespace LinearAlgebra { + template + Vector &Vector::operator= (const Vector &in_vector) + { + if (PointerComparison::equal(this, &in_vector)) + return *this; + + this->thread_loop_partitioner = in_vector.thread_loop_partitioner; + if (this->size() != in_vector.size()) + this->reinit(in_vector, true); + + dealii::internal::Vector_copy copier; + copier.dst = this->val; + copier.src = in_vector.val; + internal::parallel_for(copier, this->size(), this->thread_loop_partitioner); + + return *this; + } + + + + template + template + Vector &Vector::operator= (const Vector &in_vector) + { + this->thread_loop_partitioner = in_vector.thread_loop_partitioner; + if (this->size() != in_vector.size()) + this->reinit(in_vector, true); + + dealii::internal::Vector_copy copier; + copier.dst = this->val; + copier.src = in_vector.val; + internal::parallel_for(copier, this->size(), this->thread_loop_partitioner); + + return *this; + } + + + + template + Vector &Vector::operator= (const Number s) + { + Assert(s==static_cast(0), ExcMessage("Only 0 can be assigned to a vector.")); + (void) s; + + internal::Vector_set setter; + setter.dst = this->val; + setter.value = Number(); + + internal::parallel_for(setter, this->size(), this->thread_loop_partitioner); + + return *this; + } + + + template Vector &Vector::operator*= (const Number factor) { AssertIsFinite(factor); - const size_type size = this->size(); - for (size_type i=0; ival[i] *= factor; + + internal::Vectorization_multiply_factor vector_multiply; + vector_multiply.val = this->val; + vector_multiply.factor = factor; + + internal::parallel_for(vector_multiply, this->size(), this->thread_loop_partitioner); return *this; } @@ -59,12 +118,11 @@ namespace LinearAlgebra const Vector &down_V = dynamic_cast&>(V); Assert(down_V.size()==this->size(), ExcMessage("Cannot add two vectors with different numbers of elements")); - const size_type size = this->size(); - for (size_type i=0; ival[i] += down_V[i]; - } + + internal::Vectorization_add_v vector_add; + vector_add.val = this->val; + vector_add.v_val = down_V.val; + internal::parallel_for(vector_add, this->size(), this->thread_loop_partitioner); return *this; } @@ -81,12 +139,10 @@ namespace LinearAlgebra const Vector &down_V = dynamic_cast&>(V); Assert(down_V.size()==this->size(), ExcMessage("Cannot subtract two vectors with different numbers of elements")); - const size_type size = this->size(); - for (size_type i=0; ival[i] -= down_V[i]; - } + internal::Vectorization_subtract_v vector_subtract; + vector_subtract.val = this->val; + vector_subtract.v_val = down_V.val; + internal::parallel_for(vector_subtract, this->size(), this->thread_loop_partitioner); return *this; } @@ -105,15 +161,13 @@ namespace LinearAlgebra Assert(down_V.size()==this->size(), ExcMessage("Cannot compute the scalar product " "of two vectors with different numbers of elements")); - const size_type size = this->size(); - Number value = 0.; - for (size_type i=0; ival[i]*down_V[i]; - } + Number sum; + internal::Dot dot; + dot.X = this->val; + dot.Y = down_V.val; + internal::parallel_reduce(dot, this->size(), sum, this->thread_loop_partitioner); - return value; + return sum; } @@ -128,6 +182,20 @@ namespace LinearAlgebra + template + inline + void Vector::add(const Number a) + { + AssertIsFinite(a); + + internal::Vectorization_add_factor vector_add; + vector_add.val = this->val; + vector_add.factor = a; + internal::parallel_for(vector_add, this->size(), this->thread_loop_partitioner); + } + + + template void Vector::add(const Number a, const VectorSpaceVector &V) { @@ -140,12 +208,12 @@ namespace LinearAlgebra AssertIsFinite(a); Assert(down_V.size()==this->size(), ExcMessage("Cannot add two vectors with different numbers of elements")); - const size_type size = this->size(); - for (size_type i=0; ival[i] += a*down_V[i]; - } + + internal::Vectorization_add_av vector_add_av; + vector_add_av.val = this->val; + vector_add_av.v_val = down_V.val; + vector_add_av.factor = a; + internal::parallel_for(vector_add_av, this->size(), this->thread_loop_partitioner); } @@ -171,13 +239,14 @@ namespace LinearAlgebra AssertIsFinite(b); Assert(down_W.size()==this->size(), ExcMessage("Cannot add two vectors with different numbers of elements")); - const size_type size = this->size(); - for (size_type i=0; ival[i] += a*down_V[i]+b*down_W[i]; - } + + internal::Vectorization_add_avpbw vector_add; + vector_add.val = this->val; + vector_add.v_val = down_V.val; + vector_add.w_val = down_W.val; + vector_add.a = a; + vector_add.b = b; + internal::parallel_for(vector_add, this->size(), this->thread_loop_partitioner); } @@ -186,16 +255,21 @@ namespace LinearAlgebra void Vector::sadd(const Number s, const Number a, const VectorSpaceVector &V) { + AssertIsFinite(s); + AssertIsFinite(a); + // Check that casting will work. Assert(dynamic_cast*>(&V)!=NULL, ExcVectorTypeNotCompatible()); - *this *= s; // Downcast V. It fails, throws an exception. const Vector &down_V = dynamic_cast&>(V); - Vector tmp(down_V); - tmp *= a; - *this += tmp; + internal::Vectorization_sadd_xav vector_sadd_xav; + vector_sadd_xav.val = this->val; + vector_sadd_xav.v_val = down_V.val; + vector_sadd_xav.a = a; + vector_sadd_xav.x = s; + internal::parallel_for(vector_sadd_xav, this->size(), this->thread_loop_partitioner); } @@ -212,12 +286,11 @@ namespace LinearAlgebra dynamic_cast&>(scaling_factors); Assert(down_scaling_factors.size()==this->size(), ExcMessage("Cannot add two vectors with different numbers of elements")); - const size_type size = this->size(); - for (size_type i=0; ival[i] *= down_scaling_factors[i]; - } + + internal::Vectorization_scale vector_scale; + vector_scale.val = this->val; + vector_scale.v_val = down_scaling_factors.val; + internal::parallel_for(vector_scale, this->size(), this->thread_loop_partitioner); } @@ -225,14 +298,19 @@ namespace LinearAlgebra template void Vector::equ(const Number a, const VectorSpaceVector &V) { + AssertIsFinite(a); + // Check that casting will work. Assert(dynamic_cast*>(&V)!=NULL, ExcVectorTypeNotCompatible()); // Downcast V. If fails, throws an exception. const Vector &down_V = dynamic_cast&>(V); - *this = down_V; - *this *= a; + internal::Vectorization_equ_au vector_equ; + vector_equ.val = this->val; + vector_equ.u_val = down_V.val; + vector_equ.a = a; + internal::parallel_for(vector_equ, this->size(), this->thread_loop_partitioner); } @@ -242,10 +320,13 @@ namespace LinearAlgebra { Assert (this->size(), ExcEmptyObject()); - typename ReadWriteVector::real_type norm = 0.; - norm = l1_norm_recursive(0,this->size()-1); + typedef typename VectorSpaceVector::real_type real_type; + real_type sum; + internal::Norm1 norm1; + norm1.X = this->val; + internal::parallel_reduce(norm1, this->size(), sum, this->thread_loop_partitioner); - return norm; + return sum; } @@ -261,9 +342,12 @@ namespace LinearAlgebra // might still be finite. In that case, recompute ig (this is a rare case, // so working on the vector twice is uncritical and paid off by the extended // precision) using the BLAS approach with a weight, see e.g. dnrm2.f. - typedef typename ReadWriteVector::real_type real_type; - real_type norm_square = 0.; - norm_square = l2_norm_squared_recursive(0,this->size()-1); + typedef typename VectorSpaceVector::real_type real_type; + real_type norm_square; + internal::Norm2 norm2; + norm2.X = this->val; + internal::parallel_reduce(norm2, this->size(), norm_square, + this->thread_loop_partitioner); if (numbers::is_finite(norm_square) && norm_square>=std::numeric_limits::min()) return std::sqrt(norm_square); @@ -312,51 +396,33 @@ namespace LinearAlgebra const VectorSpaceVector &V, const VectorSpaceVector &W) { - this->add(a,V); - - return *this*W; - } - - - - template - typename VectorSpaceVector::real_type Vector::l1_norm_recursive( - unsigned int i, - unsigned int j) const - { - Assert(j>=i, ExcInternalError()); - typename ReadWriteVector::real_type norm = 0.; - - if ((j-i)!=0) - { - norm += l1_norm_recursive(i,(i+j)/2); - norm += l1_norm_recursive((i+j)/2+1,j); - } - else - norm += std::abs(this->val[i]); - - return norm; - } - - + // Check that casting will work. + Assert(dynamic_cast*>(&V)!=NULL, + ExcVectorTypeNotCompatible()); + // Check that casting will work. + Assert(dynamic_cast*>(&W)!=NULL, + ExcVectorTypeNotCompatible()); - template - typename VectorSpaceVector::real_type Vector::l2_norm_squared_recursive( - unsigned int i, - unsigned int j) const - { - Assert(j>=i, ExcInternalError()); + // Downcast V. If fails, throws an exception. + const Vector &down_V = dynamic_cast&>(V); + // Downcast W. If fails, throws an exception. + const Vector &down_W = dynamic_cast&>(W); + AssertIsFinite(a); + Assert(down_V.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); + Assert(down_W.size()==this->size(), + ExcMessage("Cannot add two vectors with different numbers of elements")); - typename ReadWriteVector::real_type norm = 0.; - if ((j-i)!=0) - { - norm += l2_norm_squared_recursive(i,(i+j)/2); - norm += l2_norm_squared_recursive((i+j)/2+1,j); - } - else - norm += std::pow(std::abs(this->val[i]),2); + Number sum; + internal::AddAndDot adder; + adder.X = this->val; + adder.a = a; + adder.V = down_V.val; + adder.W = down_W.val; + internal::parallel_reduce(adder, this->size(), sum, this->thread_loop_partitioner); + AssertIsFinite(sum); - return norm; + return sum; } diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index 695f380400..f772c2338b 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -222,6 +222,12 @@ namespace LinearAlgebra ReadWriteVector & operator= (const ReadWriteVector &in_vector); + /** + * Sets all elements of the vector to the scalar @p s. This operation is + * only allowed if @p s is equal to zero. + */ + ReadWriteVector &operator = (const Number s); + #ifdef DEAL_II_WITH_PETSC /** * Imports all the elements present in the vector's IndexSet from the input @@ -265,12 +271,6 @@ namespace LinearAlgebra std_cxx11::shared_ptr ()); #endif - /** - * Sets all elements of the vector to the scalar @p s. This operation is - * only allowed if @p s is equal to zero. - */ - ReadWriteVector &operator = (const Number s); - /** * The value returned by this function denotes the dimension of the vector * spaces that are modeled by objects of this kind. However, objects of @@ -505,8 +505,18 @@ namespace LinearAlgebra /** * Pointer to the array of local elements of this vector. */ - // TODO: use AlignedVector here for storage Number *val; + + /** + * For parallel loops with TBB, this member variable stores the affinity + * information of loops. + */ + mutable std_cxx11::shared_ptr thread_loop_partitioner; + + /** + * Make all other ReadWriteVector types friends. + */ + template friend class ReadWriteVector; }; /*@}*/ @@ -568,35 +578,6 @@ namespace LinearAlgebra - template - inline - ReadWriteVector & - ReadWriteVector::operator= (const ReadWriteVector &in_vector) - { - resize_val(in_vector.n_elements()); - stored_elements = in_vector.get_stored_elements(); - std::copy(in_vector.begin(),in_vector.end(),begin()); - - return *this; - } - - - - template - template - inline - ReadWriteVector & - ReadWriteVector::operator= (const ReadWriteVector &in_vector) - { - resize_val(in_vector.n_elements()); - stored_elements = in_vector.get_stored_elements(); - std::copy(in_vector.begin(),in_vector.end(),begin()); - - return *this; - } - - - template inline typename ReadWriteVector::size_type @@ -760,21 +741,6 @@ namespace LinearAlgebra - template - inline - ReadWriteVector & - ReadWriteVector::operator= (const Number s) - { - Assert(s==static_cast(0), ExcMessage("Only 0 can be assigned to a vector.")); - (void)s; - - std::fill(begin(),end(),Number()); - - return *this; - } - - - template template inline @@ -795,7 +761,8 @@ namespace LinearAlgebra ReadWriteVector::add (const std::vector &indices, const ReadWriteVector &values) { - for (size_type i=0; i #include -#include -#include -#include +#include + +#ifdef DEAL_II_PETSC +# include +#endif #ifdef DEAL_II_WITH_TRILINOS +# include +# include # include # include "Epetra_Import.h" #endif @@ -42,6 +46,7 @@ namespace LinearAlgebra if (val != NULL) free(val); val = NULL; + thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner()); } else { @@ -49,6 +54,8 @@ namespace LinearAlgebra free(val); Utilities::System::posix_memalign ((void **)&val, 64, sizeof(Number)*new_alloc_size); + if (new_alloc_size >= 4*internal::Vector::minimum_parallel_grain_size) + thread_loop_partitioner.reset(new parallel::internal::TBBPartitioner()); } } @@ -117,6 +124,64 @@ namespace LinearAlgebra + template + ReadWriteVector & + ReadWriteVector::operator= (const ReadWriteVector &in_vector) + { + if (PointerComparison::equal(this, &in_vector)) + return *this; + + thread_loop_partitioner = in_vector.thread_loop_partitioner; + if (n_elements() != in_vector.n_elements()) + reinit(in_vector, true); + + dealii::internal::Vector_copy copier; + copier.dst = val; + copier.src = in_vector.val; + internal::parallel_for(copier, n_elements(), thread_loop_partitioner); + + return *this; + } + + + + template + template + ReadWriteVector & + ReadWriteVector::operator= (const ReadWriteVector &in_vector) + { + thread_loop_partitioner = in_vector.thread_loop_partitioner; + if (n_elements() != in_vector.n_elements()) + reinit(in_vector, true); + + dealii::internal::Vector_copy copier; + copier.dst = val; + copier.src = in_vector.val; + internal::parallel_for(copier, n_elements(), thread_loop_partitioner); + + return *this; + } + + + + template + ReadWriteVector & + ReadWriteVector::operator= (const Number s) + { + Assert(s==static_cast(0), ExcMessage("Only 0 can be assigned to a vector.")); + (void)s; + + internal::Vector_set setter; + setter.dst = val; + setter.value = Number(); + + internal::parallel_for(setter, n_elements(), thread_loop_partitioner); + + return *this; + } + + + #ifdef DEAL_II_WITH_PETSC namespace internal { @@ -221,14 +286,16 @@ namespace LinearAlgebra { target_vector.Import(multivector, import, Insert); double *values = target_vector.Values(); - for (int i=0; i #include -#include -#include -#include -#include #include #include #include @@ -468,15 +464,6 @@ Vector::sadd (const Number x, -namespace internal -{ - namespace Vector - { - } -} - - - template template Number Vector::operator * (const Vector &v) const diff --git a/include/deal.II/lac/vector_internal.h b/include/deal.II/lac/vector_internal.h index eeedf23be4..804c84fe16 100644 --- a/include/deal.II/lac/vector_internal.h +++ b/include/deal.II/lac/vector_internal.h @@ -17,6 +17,11 @@ #ifndef dealii__vector_internal_h #define dealii__vector_internal_h +#include +#include +#include +#include + DEAL_II_NAMESPACE_OPEN namespace internal @@ -208,7 +213,7 @@ namespace internal else { DEAL_II_OPENMP_SIMD_PRAGMA - for (typename dealii::Vector::size_type i=begin; i& Vector::operator= (const Vector &) + +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + TEMPL_COPY_CONSTRUCTOR(double,float); + TEMPL_COPY_CONSTRUCTOR(float,double); + + + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + +#endif + +#undef TEMPL_COPY_CONSTRUCTOR } DEAL_II_NAMESPACE_CLOSE diff --git a/source/lac/read_write_vector.cc b/source/lac/read_write_vector.cc index d92a39f040..837b547036 100644 --- a/source/lac/read_write_vector.cc +++ b/source/lac/read_write_vector.cc @@ -20,4 +20,27 @@ DEAL_II_NAMESPACE_OPEN #include "read_write_vector.inst" +namespace LinearAlgebra +{ +// do a few functions that currently don't fit the scheme because they have +// two template arguments that need to be different (the case of same +// arguments is covered by the default copy constructor and copy operator that +// is declared separately) + +#define TEMPL_COPY_CONSTRUCTOR(S1,S2) \ + template ReadWriteVector& ReadWriteVector::operator= (const ReadWriteVector &) + +#ifndef DEAL_II_EXPLICIT_CONSTRUCTOR_BUG + TEMPL_COPY_CONSTRUCTOR(double,float); + TEMPL_COPY_CONSTRUCTOR(float,double); + + + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + TEMPL_COPY_CONSTRUCTOR(std::complex,std::complex); + +#endif + +#undef TEMPL_COPY_CONSTRUCTOR +} + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/readwritevector_assignment.cc b/tests/lac/readwritevector_assignment.cc deleted file mode 100644 index 6b358186fe..0000000000 --- a/tests/lac/readwritevector_assignment.cc +++ /dev/null @@ -1,67 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2015 by the deal.II authors -// -// This file is part of the deal.II library. -// -// The deal.II library is free software; you can use it, redistribute -// it, and/or modify it under the terms of the GNU Lesser General -// Public License as published by the Free Software Foundation; either -// version 2.1 of the License, or (at your option) any later version. -// The full text of the license can be found in the file LICENSE at -// the top level of the deal.II distribution. -// -// --------------------------------------------------------------------- - -// Check constructor and operator= - -#include "../tests.h" -#include -#include -#include - -#include - - -void test() -{ - unsigned int double_size = 2; - unsigned int float_size = 10; - IndexSet is(50); - is.add_range(0,2); - is.add_index(46); - is.add_range(10,15); - LinearAlgebra::ReadWriteVector double_vector(is); - LinearAlgebra::ReadWriteVector float_vector(float_size); - deallog << "double_size " << double_vector.n_elements() < double_vector2(double_vector); - double_vector2.print(deallog.get_file_stream()); - - for (unsigned int i=0; i