From: Bruno Turcksin Date: Tue, 24 Nov 2015 22:45:10 +0000 (-0600) Subject: Make some functions const and use covariant return types. Add ComputationPattern X-Git-Tag: v8.5.0-rc1~1180^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=531baef3cee93649d9ad67f5bd9f5f9c86231317;p=dealii.git Make some functions const and use covariant return types. Add ComputationPattern class. --- diff --git a/include/deal.II/lac/communication_pattern_base.h b/include/deal.II/lac/communication_pattern_base.h new file mode 100644 index 0000000000..9ef2852291 --- /dev/null +++ b/include/deal.II/lac/communication_pattern_base.h @@ -0,0 +1,64 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__communication_pattern_base_h +#define dealii__communication_pattern_base_h + +#include + +#ifdef DEAL_II_WITH_MPI + +#include + +DEAL_II_NAMESPACE_OPEN + +class IndexSet; + +namespace LinearAlgebra +{ + /** + * CommunicationPattern is an abstract class that is used to define a + * communication plan that can be called repeatedly to efficiently obtain + * off-processor elements. + * + * @author Bruno Turcksin, 2015. + */ + class CommunicationPatternBase + { + public: + /** + * Reinitialize the communication pattern. The first argument @p + * vector_space_vector_index_set is the index set associated to a + * VectorSpaceVector object. The second argument @p + * read_write_vector_index_set is the index set associated to a + * ReadWriteVector object. + */ + virtual void reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) = 0; + + /** + * Return a constant reference to the underlying mpi communicator. + */ + virtual const MPI_Comm &get_mpi_communicator() const = 0; + }; + +} // end of namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/include/deal.II/lac/la_vector.h b/include/deal.II/lac/la_vector.h index 09c70394af..2ce6a28c82 100644 --- a/include/deal.II/lac/la_vector.h +++ b/include/deal.II/lac/la_vector.h @@ -110,81 +110,81 @@ namespace LinearAlgebra /** * Multiply the entire vector by a fixed factor. */ - virtual VectorSpaceVector &operator*= (const Number factor); + virtual Vector &operator*= (const Number factor) override; /** * Divide the entire vector by a fixed factor. */ - virtual VectorSpaceVector &operator/= (const Number factor); + virtual Vector &operator/= (const Number factor) override; /** * Add the vector @p V to the present one. */ - virtual VectorSpaceVector &operator+= (const VectorSpaceVector &V); + virtual Vector &operator+= (const VectorSpaceVector &V) override; /** * Substract the vector @p V from the present one. */ - virtual VectorSpaceVector &operator-= (const VectorSpaceVector &V); + virtual Vector &operator-= (const VectorSpaceVector &V) override; /** * Return the scalar product of two vectors. */ - virtual Number operator* (const VectorSpaceVector &V); + virtual Number operator* (const VectorSpaceVector &V) const override; /** * Add @p a to all components. Note that @p a is a scalar not a vector. */ - virtual void add(const Number a); + virtual void add(const Number a) override; /** * Simple addition of a multiple of a vector, i.e. *this += a*V. */ - virtual void add(const Number a, const VectorSpaceVector &V); + virtual void add(const Number a, const VectorSpaceVector &V) override; /** * Multiple addition of a multiple of a vector, i.e. *this += * a*V+b*W. */ virtual void add(const Number a, const VectorSpaceVector &V, - const Number b, const VectorSpaceVector &W); + const Number b, const VectorSpaceVector &W) override; /** * Scaling and simple addition of a multiple of a vector, i.e. *this = * s*(*this)+a*V. */ virtual void sadd(const Number s, const Number a, - const VectorSpaceVector &V); + const VectorSpaceVector &V) override; /** * Scale each element of this vector by the corresponding element in the * argument. This function is mostly meant to simulate multiplication (and * immediate re-assignement) by a diagonal scaling matrix. */ - virtual void scale(const VectorSpaceVector &scaling_factors); + virtual void scale(const VectorSpaceVector &scaling_factors) override; /** * Assignement *this = a*V. */ - virtual void equ(const Number a, const VectorSpaceVector &V); + virtual void equ(const Number a, const VectorSpaceVector &V) override; /** * Return the l1 norm of the vector (i.e., the sum of the * absolute values of all entries). */ - virtual typename VectorSpaceVector::real_type l1_norm(); + virtual typename VectorSpaceVector::real_type l1_norm() const override; /** * Return the l2 norm of the vector (i.e., the square root of * the sum of the square of all entries among all processors). */ - virtual typename VectorSpaceVector::real_type l2_norm(); + virtual typename VectorSpaceVector::real_type l2_norm() const override; /** * Return the maximum norm of the vector (i.e., the maximum absolute value * among all entries and among all processors). */ - virtual typename VectorSpaceVector::real_type linfty_norm(); + virtual typename VectorSpaceVector::real_type linfty_norm() const override; /** * Perform a combined operation of a vector addition and a subsequent @@ -197,13 +197,13 @@ namespace LinearAlgebra */ virtual Number add_and_dot(const Number a, const VectorSpaceVector &V, - const VectorSpaceVector &W); + const VectorSpaceVector &W) override; /** * Return the global size of the vector, equal to the sum of the number of * locally owned indices among all processors. */ - size_type size() const; + virtual size_type size() const override; /** * Return an index set that describes which elements of this vector are @@ -216,7 +216,7 @@ namespace LinearAlgebra * vec.locally_owned_elements() == complete_index_set(vec.size()) * @endcode */ - const dealii::IndexSet locally_owned_elements() const; + virtual dealii::IndexSet locally_owned_elements() const override; /** * Prints the vector to the output stream @p out. @@ -224,7 +224,7 @@ namespace LinearAlgebra virtual void print(std::ostream &out, const unsigned int precision=3, const bool scientific=true, - const bool across=true) const; + const bool across=true) const override; /** * Write the vector en bloc to a file. This is done in a binary mode, so * the output is neither readable by humans nor (probably) by other @@ -248,7 +248,7 @@ namespace LinearAlgebra /** * Returns the memory consumption of this class in bytes. */ - std::size_t memory_consumption() const; + virtual std::size_t memory_consumption() const override; /** * Attempt to perform an operation between two incompatible vector types. @@ -264,7 +264,7 @@ namespace LinearAlgebra * large vector. */ typename VectorSpaceVector::real_type l1_norm_recursive(unsigned int i, - unsigned int j); + unsigned int j) const; /** * Compute the squared L2 norm in a recursive way by dividing the vector @@ -273,7 +273,7 @@ namespace LinearAlgebra */ typename VectorSpaceVector::real_type l2_norm_squared_recursive( unsigned int i, - unsigned int j); + unsigned int j) const; /** * Serialize the data of this object using boost. This function is @@ -391,9 +391,9 @@ namespace LinearAlgebra template inline - const dealii::IndexSet Vector::locally_owned_elements() const + dealii::IndexSet Vector::locally_owned_elements() const { - return ReadWriteVector::get_stored_elements(); + return IndexSet(ReadWriteVector::get_stored_elements()); } diff --git a/include/deal.II/lac/la_vector.templates.h b/include/deal.II/lac/la_vector.templates.h index 95e8465435..bc014e9dae 100644 --- a/include/deal.II/lac/la_vector.templates.h +++ b/include/deal.II/lac/la_vector.templates.h @@ -25,7 +25,7 @@ DEAL_II_NAMESPACE_OPEN namespace LinearAlgebra { template - VectorSpaceVector &Vector::operator*= (const Number factor) + Vector &Vector::operator*= (const Number factor) { AssertIsFinite(factor); for (unsigned int i=0; isize(); ++i) @@ -37,7 +37,7 @@ namespace LinearAlgebra template - VectorSpaceVector &Vector::operator/= (const Number factor) + Vector &Vector::operator/= (const Number factor) { AssertIsFinite(factor); Assert(factor!=Number(0.), ExcZero()); @@ -49,7 +49,7 @@ namespace LinearAlgebra template - VectorSpaceVector &Vector::operator+= (const VectorSpaceVector &V) + Vector &Vector::operator+= (const VectorSpaceVector &V) { // Check that casting will work. Assert(dynamic_cast*>(&V)!=NULL, ExcVectorTypeNotCompatible()); @@ -70,7 +70,7 @@ namespace LinearAlgebra template - VectorSpaceVector &Vector::operator-= (const VectorSpaceVector &V) + Vector &Vector::operator-= (const VectorSpaceVector &V) { // Check that casting will work. Assert(dynamic_cast*>(&V)!=NULL, ExcVectorTypeNotCompatible()); @@ -91,7 +91,7 @@ namespace LinearAlgebra template - Number Vector::operator* (const VectorSpaceVector &V) + Number Vector::operator* (const VectorSpaceVector &V) const { // Check that casting will work. Assert(dynamic_cast*>(&V)!=NULL, @@ -221,7 +221,7 @@ namespace LinearAlgebra template - typename VectorSpaceVector::real_type Vector::l1_norm() + typename VectorSpaceVector::real_type Vector::l1_norm() const { Assert (this->size(), ExcEmptyObject()); @@ -235,7 +235,7 @@ namespace LinearAlgebra template - typename VectorSpaceVector::real_type Vector::l2_norm() + typename VectorSpaceVector::real_type Vector::l2_norm() const { Assert (this->size(), ExcEmptyObject()); @@ -277,7 +277,7 @@ namespace LinearAlgebra template - typename VectorSpaceVector::real_type Vector::linfty_norm() + typename VectorSpaceVector::real_type Vector::linfty_norm() const { typename ReadWriteVector::real_type norm = 0.; for (unsigned int i=0; isize(); ++i) @@ -303,7 +303,7 @@ namespace LinearAlgebra template typename VectorSpaceVector::real_type Vector::l1_norm_recursive( unsigned int i, - unsigned int j) + unsigned int j) const { Assert(j>=i, ExcInternalError()); typename ReadWriteVector::real_type norm = 0.; @@ -324,7 +324,7 @@ namespace LinearAlgebra template typename VectorSpaceVector::real_type Vector::l2_norm_squared_recursive( unsigned int i, - unsigned int j) + unsigned int j) const { Assert(j>=i, ExcInternalError()); diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index a2574caa07..f39e50c8de 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -35,6 +35,11 @@ DEAL_II_NAMESPACE_OPEN +namespace LinearAlgebra +{ + class CommunicationPatternBase; +} + #ifdef DEAL_II_WITH_PETSC namespace PETScWrappers { @@ -218,27 +223,42 @@ namespace LinearAlgebra #ifdef DEAL_II_WITH_PETSC /** - * Imports all the elements present in the vector's IndexSet from the - * input vector @p petsc_vec. - */ - ReadWriteVector & - operator= (const PETScWrappers::MPI::Vector &petsc_vec); + * Imports all the elements present in the vector's IndexSet from the input + * vector @p petsc_vec. VectorOperation::values @p operation is used to decide + * if the elements in @p V should be added to the current vector or replace + * the current elements. The last parameter can be used if the same + * communication pattern is used multiple times. This can be used to improve + * performance. + */ + void import(const PETScWrappers::MPI::Vector &petsc_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL); #endif #ifdef DEAL_II_WITH_TRILINOS - /** - * Imports all the elements present in the vector's IndexSet from the - * input vector @p trilinos_vec. - */ - ReadWriteVector & - operator= (const TrilinosWrappers::MPI::Vector &trilinos_vec); - /** * Imports all the elements present in the vector's IndexSet from the input - * vector @p epetra_vec. - */ - ReadWriteVector & - operator= (const EpetraWrappers::Vector &epetra_vec); + * vector @p trilinos_vec. VectorOperation::values @p operation is used to + * decide if the elements in @p V should be added to the current vector or + * replace the current elements. The last parameter can be used if the same + * communication pattern is used multiple times. This can be used to improve + * performance. + */ + void import(const TrilinosWrappers::MPI::Vector &trilinos_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL); + + /** + * imports all the elements present in the vector's IndexSet from the input + * vector @p epetra_vec. VectorOperation::values @p operation is used to + * decide if the elements in @p V should be added to the current vector or + * replace the current elements. The last parameter can be used if the same + * communication pattern is used multiple times. This can be used to improve + * performance. + */ + void import(const EpetraWrappers::Vector &epetra_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL); #endif /** @@ -426,10 +446,14 @@ namespace LinearAlgebra #ifdef DEAL_II_WITH_TRILINOS /** * Imports all the elements present in the vector's IndexSet from the input - * vector @p multivector. - */ - ReadWriteVector &import(const Epetra_MultiVector &multivector, - const IndexSet &locally_owned_elements); + * vector @p multivector. This is an helper function and it should not be + * used directly. + */ + void import(const Epetra_MultiVector &multivector, + const IndexSet &locally_owned_elements, + VectorOperation::values operation, + const MPI_Comm &mpi_comm, + const CommunicationPatternBase *communication_pattern); #endif /** @@ -448,11 +472,27 @@ namespace LinearAlgebra */ void resize_val (const size_type new_allocated_size); + /** + * + */ + void create_epetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm); + /** * Indices of the elements stored. */ IndexSet stored_elements; + /** + * Indices of the elements stored on a VectorSpaceVector + */ + std::shared_ptr source_stored_elements; + + /** + * + */ + std::shared_ptr comm_pattern; + /** * Pointer to the array of local elements of this vector. */ @@ -471,7 +511,9 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector () : - val (NULL) + source_stored_elements(nullptr), + comm_pattern(nullptr), + val(nullptr) {} @@ -481,7 +523,9 @@ namespace LinearAlgebra ReadWriteVector::ReadWriteVector (const ReadWriteVector &v) : Subscriptor(), - val (NULL) + source_stored_elements(nullptr), + comm_pattern(nullptr), + val(nullptr) { this->operator=(v); } @@ -492,7 +536,9 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector (const size_type size) : - val (NULL) + source_stored_elements(nullptr), + comm_pattern(nullptr), + val(nullptr) { reinit (size, false); } @@ -503,7 +549,9 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector (const IndexSet &locally_stored_indices) : - val (NULL) + source_stored_elements(nullptr), + comm_pattern(nullptr), + val(nullptr) { reinit (locally_stored_indices); } diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index 936574e918..0e42ecf04b 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -1,6 +1,6 @@ // --------------------------------------------------------------------- // -// Copyright (C) 2015, 2016 by the deal.II authors +// Copyright (C) 2015-2016 by the deal.II authors // // This file is part of the deal.II library. // @@ -24,7 +24,8 @@ #include #ifdef DEAL_II_WITH_TRILINOS -# include +# include +# include "Epetra_Import.h" #endif DEAL_II_NAMESPACE_OPEN @@ -38,13 +39,13 @@ namespace LinearAlgebra { if (new_alloc_size == 0) { - if (val != NULL) + if (val != nullptr) free(val); - val = NULL; + val = nullptr; } else { - if (val != NULL) + if (val != nullptr) free(val); Utilities::System::posix_memalign ((void **)&val, 64, sizeof(Number)*new_alloc_size); @@ -67,6 +68,10 @@ namespace LinearAlgebra // set entries to zero if so requested if (omit_zeroing_entries == false) this->operator = (Number()); + + // reset the communication patter + source_stored_elements = nullptr; + comm_pattern = nullptr; } @@ -83,6 +88,10 @@ namespace LinearAlgebra if (omit_zeroing_entries == false) this->operator= (Number()); + + // reset the communication patter + source_stored_elements = nullptr; + comm_pattern = nullptr; } @@ -100,6 +109,10 @@ namespace LinearAlgebra // initialize to zero if (omit_zeroing_entries == false) this->operator= (Number()); + + // reset the communication patter + source_stored_elements = nullptr; + comm_pattern = nullptr; } @@ -135,8 +148,10 @@ namespace LinearAlgebra template - ReadWriteVector & - ReadWriteVector::operator = (const PETScWrappers::MPI::Vector &petsc_vec) + void + ReadWriteVector::import(const PETScWrappers::MPI::Vector &petsc_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern) { //TODO: this works only if no communication is needed. Assert(petsc_vec.locally_owned_elements() == stored_elements, @@ -153,10 +168,6 @@ namespace LinearAlgebra // restore the representation of the vector ierr = VecRestoreArray (static_cast(petsc_vec), &start_ptr); AssertThrow (ierr == 0, ExcPETScError(ierr)); - - // return a pointer to this object per normal c++ operator overloading - // semantics - return *this; } #endif @@ -164,83 +175,83 @@ namespace LinearAlgebra #ifdef DEAL_II_WITH_TRILINOS template - ReadWriteVector & - ReadWriteVector::import(const Epetra_MultiVector &multivector, - const IndexSet &locally_owned_elements) + void + ReadWriteVector::import(const Epetra_MultiVector &multivector, + const IndexSet &source_elements, + VectorOperation::values operation, + const MPI_Comm &mpi_comm, + const CommunicationPatternBase *communication_pattern) { - // Copy the local elements - std::vector trilinos_indices, readwrite_indices; - locally_owned_elements.fill_index_vector(trilinos_indices); - stored_elements.fill_index_vector(readwrite_indices); - std::vector intersection(trilinos_indices.size()); - std::vector::iterator end; - end = std::set_intersection(trilinos_indices.begin(), trilinos_indices.end(), - readwrite_indices.begin(), readwrite_indices.end(), - intersection.begin()); - const size_t intersection_size = end-intersection.begin(); - intersection.resize(intersection_size); - std::vector trilinos_position(intersection_size); - unsigned int j=0; - for (size_t i=0; i (comm_pattern.get()); + if (epetra_comm_pattern == nullptr) + create_epetra_comm_pattern(source_elements, mpi_comm); + } + else + create_epetra_comm_pattern(source_elements, mpi_comm); + + epetra_comm_pattern = + dynamic_cast (comm_pattern.get()); + } + else + { + epetra_comm_pattern = + dynamic_cast (communication_pattern); + AssertThrow(epetra_comm_pattern != nullptr, + ExcMessage(std::string("The communication pattern is not of type ") + + "LinearAlgebra::EpetraWrappers::CommunicationPattern.")); } - double *values = multivector.Values(); - for (size_t i=0; iget_epetra_import()); -#ifdef DEAL_II_WITH_MPI - // Copy the off-processor elements if necessary - if (intersection_size != n_elements()) + Epetra_FEVector target_vector(import.TargetMap()); + + if (operation==VectorOperation::insert) { - Epetra_BlockMap source_map = multivector.Map(); - std::vector off_proc_indices(readwrite_indices.size()); - // Subtract the elements of readwrite that are in intersection. - end = std::set_difference(readwrite_indices.begin(), readwrite_indices.end(), - intersection.begin(), intersection.end(), off_proc_indices.begin()); - off_proc_indices.resize(end-off_proc_indices.begin()); - // Cast size_type to TrilinosWrappers::type::int_type - TrilinosWrappers::types::int_type *off_proc_array = - new TrilinosWrappers::types::int_type [off_proc_indices.size()]; - for (size_t i=0; i ( - off_proc_indices[i]); - Epetra_Map target_map(off_proc_indices.size(),off_proc_indices.size(), - off_proc_array,0,source_map.Comm()); - delete [] off_proc_array; - Epetra_Import import(target_map, source_map); - Epetra_FEVector target_vector(target_map); target_vector.Import(multivector, import, Insert); - values = target_vector.Values(); - for (unsigned int i=0; i - ReadWriteVector & - ReadWriteVector::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec) + void + ReadWriteVector::import(const TrilinosWrappers::MPI::Vector &trilinos_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern) { - return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements()); + import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements(), + operation, trilinos_vec.get_mpi_communicator(), communication_pattern); } + template - ReadWriteVector & - ReadWriteVector::operator = (const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec) + void + ReadWriteVector::import(const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern) { - return import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements()); + import(trilinos_vec.trilinos_vector(), trilinos_vec.locally_owned_elements(), + operation, trilinos_vec.get_mpi_communicator(), communication_pattern); } +#endif + template @@ -296,6 +307,18 @@ namespace LinearAlgebra out.flags (old_flags); out.precision(old_precision); } + + + + template + void + ReadWriteVector::create_epetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm) + { + source_stored_elements = std::make_shared(source_index_set); + comm_pattern.reset(new EpetraWrappers::CommunicationPattern( + *source_stored_elements, stored_elements, mpi_comm)); + } } // end of namespace LinearAlgebra diff --git a/include/deal.II/lac/trilinos_epetra_communication_pattern.h b/include/deal.II/lac/trilinos_epetra_communication_pattern.h new file mode 100644 index 0000000000..d25f9187d1 --- /dev/null +++ b/include/deal.II/lac/trilinos_epetra_communication_pattern.h @@ -0,0 +1,90 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015-2016 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. +// +// --------------------------------------------------------------------- + +#ifndef dealii__trilinos_epetra_communication_pattern_h +#define dealii__trilinos_epetra_communication_pattern_h + + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_MPI + +#include +#include +#include "Epetra_Import.h" + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace EpetraWrappers + { + /** + * This class implements a wrapper to Trilinos Import. + */ + class CommunicationPattern : public CommunicationPatternBase + { + public: + /** + * Reinitialize the communication pattern. The first argument @p + * vector_space_vector_index_set is the index set associated to a + * VectorSpaceVector object. The second argument @p + * read_write_vector_index_set is the index set associated to a + * ReadWriteVector object. + */ + CommunicationPattern(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator); + + /** + * Reinitialize the object. + */ + void reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) override; + + /** + * Return the underlying MPI communicator. + */ + const MPI_Comm &get_mpi_communicator() const override; + + /** + * Return the underlying Epetra_Import object. + */ + const Epetra_Import &get_epetra_import() const; + + private: + /** + * Shared pointer to the MPI communicator used. + */ + std_cxx11::shared_ptr comm; + + /** + * Shared pointer to the Epetra_Import object used. + */ + std_cxx11::shared_ptr import; + }; + } // end of namespace EpetraWrappers +} // end of namespace LinearAlgebra + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif + +#endif diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h index 61301b6f6b..527a9af254 100644 --- a/include/deal.II/lac/trilinos_epetra_vector.h +++ b/include/deal.II/lac/trilinos_epetra_vector.h @@ -41,7 +41,7 @@ namespace LinearAlgebra namespace EpetraWrappers { /** - * This class implements a wrapper to the the Trilinos distributed vector + * This class implements a wrapper to the Trilinos distributed vector * class Epetra_FEVector. This class is derived from the * LinearAlgebra::VectorSpaceVector class. Note however that Epetra only * works with Number = double. @@ -67,11 +67,11 @@ namespace LinearAlgebra /** * This constructor takes an IndexSet that defines how to distribute the * individual components among the MPI processors. Since it also - * includes information abtou the size of the vector, this is all we + * includes information about the size of the vector, this is all we * need to genereate a %parallel vector. */ - Vector(const IndexSet ¶llel_partitioner, - const MPI_Comm &communicator = MPI_COMM_WORLD); + explicit Vector(const IndexSet ¶llel_partitioner, + const MPI_Comm &communicator); /** * Reinit functionality. This function destroys the old vector content @@ -80,7 +80,7 @@ namespace LinearAlgebra * filled with zero (false) or left untouched (true). */ void reinit (const IndexSet ¶llel_partitioner, - const MPI_Comm &communicator = MPI_COMM_WORLD, + const MPI_Comm &communicator, const bool omit_zeroing_entries = false); /** @@ -100,28 +100,28 @@ namespace LinearAlgebra /** * Multiply the entire vector by a fixed factor. */ - virtual VectorSpaceVector &operator*= (const double factor); + virtual Vector &operator*= (const double factor); /** * Divide the entire vector by a fixed factor. */ - virtual VectorSpaceVector &operator/= (const double factor); + virtual Vector &operator/= (const double factor); /** * Add the vector @p V to the present one. */ - virtual VectorSpaceVector &operator+= (const VectorSpaceVector &V); + virtual Vector &operator+= (const VectorSpaceVector &V); /** * Substract the vector @p V from the present one. */ - virtual VectorSpaceVector &operator-= (const VectorSpaceVector &V); + virtual Vector &operator-= (const VectorSpaceVector &V); /** * Return the scalar product of two vectors. The vectors need to have the * same layout. */ - virtual double operator* (const VectorSpaceVector &V); + virtual double operator* (const VectorSpaceVector &V) const; /** * Add @p a to all components. Note that @p is a scalar not a vector. @@ -166,19 +166,19 @@ namespace LinearAlgebra * Returns the l1 norm of the vector (i.e., the sum of the * absolute values of all entries among all processors). */ - virtual double l1_norm(); + virtual double l1_norm() const; /** * Returns the l2 norm of the vector (i.e., the square root of * the sum of the square of all entries among all processors). */ - virtual double l2_norm(); + virtual double l2_norm() const; /** * Returns the maximum norm of the vector (i.e., the maximum absolute value * among all entries and among all processors). */ - virtual double linfty_norm(); + virtual double linfty_norm() const; /** * Performs a combined operation of a vector addition and a subsequent @@ -221,7 +221,7 @@ namespace LinearAlgebra * vec.locally_owned_elements() == complete_index_set(vec.size()) * @endcode */ - virtual const ::dealii::IndexSet locally_owned_elements() const; + virtual ::dealii::IndexSet locally_owned_elements() const; /** * Return a const reference to the underlying Trilinos diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index ec008f89d3..26c0b584e2 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -69,10 +69,22 @@ namespace LinearAlgebra */ virtual VectorSpaceVector &operator-= (const VectorSpaceVector &V) = 0; + /** + * Imports all the elements present in the vector's IndexSet from the input + * vector @p V. VectorOperation::values @p operation is used to decide if + * the elements in @p V should be added to the current vector or replace the + * current elements. The last parameter can be used if the same + * communication pattern is used multiple times. This can be used to improve + * performance. + */ + virtual void Import(const ReadWriteVector &V, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL); + /** * Return the scalar product of two vectors. */ - virtual Number operator* (const VectorSpaceVector &V) = 0; + virtual Number operator* (const VectorSpaceVector &V) const = 0; /** * Add @p a to all components. Note that @p a is a scalar not a vector. @@ -113,19 +125,19 @@ namespace LinearAlgebra * Return the l1 norm of the vector (i.e., the sum of the * absolute values of all entries among all processors). */ - virtual real_type l1_norm() = 0; + virtual real_type l1_norm() const = 0; /** * Return the l2 norm of the vector (i.e., the square root of * the sum of the square of all entries among all processors). */ - virtual real_type l2_norm() = 0; + virtual real_type l2_norm() const = 0; /** * Return the maximum norm of the vector (i.e., the maximum absolute value * among all entries and among all processors). */ - virtual real_type linfty_norm() = 0; + virtual real_type linfty_norm() const = 0; /** * Perform a combined operation of a vector addition and a subsequent @@ -164,7 +176,7 @@ namespace LinearAlgebra * vec.locally_owned_elements() == complete_index_set(vec.size()) * @endcode */ - virtual const dealii::IndexSet locally_owned_elements() const = 0; + virtual dealii::IndexSet locally_owned_elements() const = 0; /** * Print the vector to the output stream @p out. diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt index b131f2746f..9dc145c2c5 100644 --- a/source/lac/CMakeLists.txt +++ b/source/lac/CMakeLists.txt @@ -111,6 +111,7 @@ IF(DEAL_II_WITH_TRILINOS) ${_src} trilinos_block_sparse_matrix.cc trilinos_block_vector.cc + trilinos_epetra_communication_pattern.cc trilinos_epetra_vector.cc trilinos_parallel_block_vector.cc trilinos_precondition.cc diff --git a/source/lac/trilinos_epetra_communication_pattern.cc b/source/lac/trilinos_epetra_communication_pattern.cc new file mode 100644 index 0000000000..ef93fc9db1 --- /dev/null +++ b/source/lac/trilinos_epetra_communication_pattern.cc @@ -0,0 +1,75 @@ +// --------------------------------------------------------------------- +// +// 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. +// +// --------------------------------------------------------------------- + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +#ifdef DEAL_II_WITH_MPI + +#include +#include "Epetra_Map.h" + +DEAL_II_NAMESPACE_OPEN + +namespace LinearAlgebra +{ + namespace EpetraWrappers + { + CommunicationPattern::CommunicationPattern(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) + { + reinit(vector_space_vector_index_set, read_write_vector_index_set, communicator); + } + + + + void CommunicationPattern::reinit(const IndexSet &vector_space_vector_index_set, + const IndexSet &read_write_vector_index_set, + const MPI_Comm &communicator) + { + comm = std::make_shared(communicator); + + Epetra_Map vector_space_vector_map = vector_space_vector_index_set.make_trilinos_map(*comm); + Epetra_Map read_write_vector_map = read_write_vector_index_set.make_trilinos_map(*comm); + + // Target map is read_write_vector_map + // Source map is vector_space_vector_map. This map must have uniquely + // owned GID. + import.reset(new Epetra_Import(read_write_vector_map, vector_space_vector_map)); + } + + + + const MPI_Comm &CommunicationPattern::get_mpi_communicator() const + { + return *comm; + } + + + + const Epetra_Import &CommunicationPattern::get_epetra_import() const + { + return *import; + } + } +} + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index 4bf7b32ea9..6340d56fb9 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -16,6 +16,7 @@ #include #ifdef DEAL_II_WITH_TRILINOS + #ifdef DEAL_II_WITH_MPI #include @@ -184,7 +185,7 @@ namespace LinearAlgebra - VectorSpaceVector &Vector::operator*= (const double factor) + Vector &Vector::operator*= (const double factor) { AssertIsFinite(factor); vector->Scale(factor); @@ -194,7 +195,7 @@ namespace LinearAlgebra - VectorSpaceVector &Vector::operator/= (const double factor) + Vector &Vector::operator/= (const double factor) { AssertIsFinite(factor); Assert(factor!=0., ExcZero()); @@ -205,7 +206,7 @@ namespace LinearAlgebra - VectorSpaceVector &Vector::operator+= (const VectorSpaceVector &V) + Vector &Vector::operator+= (const VectorSpaceVector &V) { // Check that casting will work. Assert(dynamic_cast(&V)!=NULL, ExcVectorTypeNotCompatible()); @@ -247,7 +248,7 @@ namespace LinearAlgebra - VectorSpaceVector &Vector::operator-= (const VectorSpaceVector &V) + Vector &Vector::operator-= (const VectorSpaceVector &V) { this->add(-1.,V); @@ -256,7 +257,7 @@ namespace LinearAlgebra - double Vector::operator* (const VectorSpaceVector &V) + double Vector::operator* (const VectorSpaceVector &V) const { // Check that casting will work. Assert(dynamic_cast(&V)!=NULL, @@ -392,7 +393,7 @@ namespace LinearAlgebra - double Vector::l1_norm() + double Vector::l1_norm() const { double norm(0.); int ierr = vector->Norm1(&norm); @@ -403,7 +404,7 @@ namespace LinearAlgebra - double Vector::l2_norm() + double Vector::l2_norm() const { double norm(0.); int ierr = vector->Norm2(&norm); @@ -414,7 +415,7 @@ namespace LinearAlgebra - double Vector::linfty_norm() + double Vector::linfty_norm() const { double norm(0.); int ierr = vector->NormInf(&norm); @@ -456,7 +457,7 @@ namespace LinearAlgebra - const ::dealii::IndexSet Vector::locally_owned_elements() const + ::dealii::IndexSet Vector::locally_owned_elements() const { const Epetra_Map *map = dynamic_cast(&(vector->Map())); return IndexSet(*map); @@ -499,10 +500,10 @@ namespace LinearAlgebra out.setf(std::ios::fixed, std::ios::floatfield); if (across) - for (size_type i=0; iMyLength(); ++i) out << val[i] << ' '; else - for (size_type i=0; iMyLength(); ++i) out << val[i] << std::endl; out << std::endl; diff --git a/tests/trilinos/epetra_vector_01.cc b/tests/trilinos/epetra_vector_01.cc index c8fe4ee422..3059910224 100644 --- a/tests/trilinos/epetra_vector_01.cc +++ b/tests/trilinos/epetra_vector_01.cc @@ -44,14 +44,14 @@ void test() parallel_partitioner_1.compress(); parallel_partitioner_2.compress(); LinearAlgebra::EpetraWrappers::Vector a; - LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1); + LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1, MPI_COMM_WORLD); LinearAlgebra::EpetraWrappers::Vector c(b); AssertThrow(a.size()==0, ExcMessage("Vector has the wrong size.")); AssertThrow(b.size()==10, ExcMessage("Vector has the wrong size.")); AssertThrow(c.size()==10, ExcMessage("Vector has the wrong size.")); - a.reinit(parallel_partitioner_2); + a.reinit(parallel_partitioner_2, MPI_COMM_WORLD); AssertThrow(a.size()==10, ExcMessage("Vector has the wrong size.")); AssertThrow(parallel_partitioner_1==b.locally_owned_elements(), diff --git a/tests/trilinos/epetra_vector_02.cc b/tests/trilinos/epetra_vector_02.cc index 746de05288..687038a393 100644 --- a/tests/trilinos/epetra_vector_02.cc +++ b/tests/trilinos/epetra_vector_02.cc @@ -40,9 +40,9 @@ void test() } parallel_partitioner_1.compress(); parallel_partitioner_2.compress(); - LinearAlgebra::EpetraWrappers::Vector a(parallel_partitioner_1); - LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1); - LinearAlgebra::EpetraWrappers::Vector c(parallel_partitioner_2); + LinearAlgebra::EpetraWrappers::Vector a(parallel_partitioner_1, MPI_COMM_WORLD); + LinearAlgebra::EpetraWrappers::Vector b(parallel_partitioner_1, MPI_COMM_WORLD); + LinearAlgebra::EpetraWrappers::Vector c(parallel_partitioner_2, MPI_COMM_WORLD); IndexSet read_write_index_set(10); if (rank==0)