From 37371b55ec27f17949c4e39fb0a1165b00e60faa Mon Sep 17 00:00:00 2001 From: Bruno Turcksin Date: Fri, 12 Feb 2016 14:13:36 -0500 Subject: [PATCH] Add import function to LinearAlgebra::TrilinosWrappers::Vector. --- .../deal.II/lac/communication_pattern_base.h | 5 +- include/deal.II/lac/read_write_vector.h | 24 ++-- .../deal.II/lac/read_write_vector.templates.h | 46 +++--- include/deal.II/lac/trilinos_epetra_vector.h | 33 ++++- include/deal.II/lac/vector_space_vector.h | 2 +- source/lac/trilinos_epetra_vector.cc | 131 +++++++----------- 6 files changed, 123 insertions(+), 118 deletions(-) diff --git a/include/deal.II/lac/communication_pattern_base.h b/include/deal.II/lac/communication_pattern_base.h index 9ef2852291..e63a2d25fd 100644 --- a/include/deal.II/lac/communication_pattern_base.h +++ b/include/deal.II/lac/communication_pattern_base.h @@ -31,7 +31,10 @@ 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. + * off-processor elements. The idea is to decouple the communication pattern + * from the vectors. The goal is to reuse the same communication pattern for + * different vectors. This is similar to the way SparseMatrix and + * SparsityPattern works. * * @author Bruno Turcksin, 2015. */ diff --git a/include/deal.II/lac/read_write_vector.h b/include/deal.II/lac/read_write_vector.h index f39e50c8de..cb61c2c9a7 100644 --- a/include/deal.II/lac/read_write_vector.h +++ b/include/deal.II/lac/read_write_vector.h @@ -30,6 +30,7 @@ #include #ifdef DEAL_II_WITH_TRILINOS +#include #include "Epetra_MultiVector.h" #endif @@ -249,7 +250,7 @@ namespace LinearAlgebra const CommunicationPatternBase *communication_pattern = NULL); /** - * imports all the elements present in the vector's IndexSet from the input + * 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 @@ -472,11 +473,15 @@ namespace LinearAlgebra */ void resize_val (const size_type new_allocated_size); +#ifdef DEAL_II_WITH_TRILINOS /** - * + * Return a EpetraWrappers::Communication pattern and store it for future + * use. */ - void create_epetra_comm_pattern(const IndexSet &source_index_set, - const MPI_Comm &mpi_comm); + EpetraWrappers::CommunicationPattern + create_epetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm); +#endif /** * Indices of the elements stored. @@ -484,12 +489,13 @@ namespace LinearAlgebra IndexSet stored_elements; /** - * Indices of the elements stored on a VectorSpaceVector + * IndexSet of the elements of the last imported vector; */ - std::shared_ptr source_stored_elements; + IndexSet source_stored_elements; /** - * + * CommunicationPattern for the communication between the + * source_stored_elements IndexSet and the current vector. */ std::shared_ptr comm_pattern; @@ -511,7 +517,6 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector () : - source_stored_elements(nullptr), comm_pattern(nullptr), val(nullptr) {} @@ -523,7 +528,6 @@ namespace LinearAlgebra ReadWriteVector::ReadWriteVector (const ReadWriteVector &v) : Subscriptor(), - source_stored_elements(nullptr), comm_pattern(nullptr), val(nullptr) { @@ -536,7 +540,6 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector (const size_type size) : - source_stored_elements(nullptr), comm_pattern(nullptr), val(nullptr) { @@ -549,7 +552,6 @@ namespace LinearAlgebra inline ReadWriteVector::ReadWriteVector (const IndexSet &locally_stored_indices) : - source_stored_elements(nullptr), comm_pattern(nullptr), val(nullptr) { diff --git a/include/deal.II/lac/read_write_vector.templates.h b/include/deal.II/lac/read_write_vector.templates.h index 0e42ecf04b..e3a7ea8a56 100644 --- a/include/deal.II/lac/read_write_vector.templates.h +++ b/include/deal.II/lac/read_write_vector.templates.h @@ -70,7 +70,7 @@ namespace LinearAlgebra this->operator = (Number()); // reset the communication patter - source_stored_elements = nullptr; + source_stored_elements.clear(); comm_pattern = nullptr; } @@ -90,7 +90,7 @@ namespace LinearAlgebra this->operator= (Number()); // reset the communication patter - source_stored_elements = nullptr; + source_stored_elements.clear(); comm_pattern = nullptr; } @@ -111,7 +111,7 @@ namespace LinearAlgebra this->operator= (Number()); // reset the communication patter - source_stored_elements = nullptr; + source_stored_elements.clear(); comm_pattern = nullptr; } @@ -182,27 +182,31 @@ namespace LinearAlgebra const MPI_Comm &mpi_comm, const CommunicationPatternBase *communication_pattern) { - const EpetraWrappers::CommunicationPattern *epetra_comm_pattern = nullptr; + std_cxx11::shared_ptr epetra_comm_pattern; + // If no communication pattern is given, create one. Otherwise, use the one + // given. if (communication_pattern == nullptr) { - if (source_elements == *source_stored_elements) + // The first time import is called, we create a communication pattern. + // Check if the communication pattern already exists and if it can be + // reused. + if (source_elements == source_stored_elements) { - epetra_comm_pattern = - dynamic_cast (comm_pattern.get()); + epetra_comm_pattern.reset( + dynamic_cast (comm_pattern.get())); if (epetra_comm_pattern == nullptr) - create_epetra_comm_pattern(source_elements, mpi_comm); + epetra_comm_pattern = std::make_shared( + 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()); + epetra_comm_pattern = std::make_shared( + create_epetra_comm_pattern(source_elements, mpi_comm)); } else { - epetra_comm_pattern = - dynamic_cast (communication_pattern); + epetra_comm_pattern.reset( + dynamic_cast (communication_pattern)); AssertThrow(epetra_comm_pattern != nullptr, ExcMessage(std::string("The communication pattern is not of type ") + "LinearAlgebra::EpetraWrappers::CommunicationPattern.")); @@ -310,15 +314,21 @@ namespace LinearAlgebra +#ifdef DEAL_II_WITH_TRILINOS template - void + EpetraWrappers::CommunicationPattern 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)); + source_stored_elements = source_index_set; + EpetraWrappers::CommunicationPattern epetra_comm_pattern( + source_stored_elements, stored_elements, mpi_comm); + comm_pattern = std::make_shared( + epetra_comm_pattern); + + return epetra_comm_pattern; } +#endif } // end of namespace LinearAlgebra diff --git a/include/deal.II/lac/trilinos_epetra_vector.h b/include/deal.II/lac/trilinos_epetra_vector.h index 527a9af254..67fd2efe71 100644 --- a/include/deal.II/lac/trilinos_epetra_vector.h +++ b/include/deal.II/lac/trilinos_epetra_vector.h @@ -27,6 +27,7 @@ #include #include #include +#include #include "mpi.h" #include "Epetra_FEVector.h" @@ -91,11 +92,17 @@ namespace LinearAlgebra Vector &operator= (const Vector &V); /** - * Copy function. Copy the elements from the ReadWriteVector @p V in the - * Vector. If elements if @p V are not locally owned, there is a - * communication. + * 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. */ - Vector &operator= (const ::dealii::LinearAlgebra::ReadWriteVector &V); + virtual void import(const ReadWriteVector &V, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern = NULL) + override; /** * Multiply the entire vector by a fixed factor. @@ -272,10 +279,28 @@ namespace LinearAlgebra << " occurred while calling a Trilinos function"); private: + /** + * Create the CommunicationPattern for the communication between the + * IndexSet @param source_index_set and the current vector. + */ + void create_epetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm); + /** * Pointer to the actual Epetra vector object. */ std_cxx11::shared_ptr vector; + + /** + * IndexSet of the elements of the last imported vector. + */ + ::dealii::IndexSet source_stored_elements; + + /** + * CommunicationPattern for the communication between the + * source_stored_elements IndexSet and the current vector. + */ + std_cxx11::shared_ptr epetra_comm_pattern; }; } } diff --git a/include/deal.II/lac/vector_space_vector.h b/include/deal.II/lac/vector_space_vector.h index 26c0b584e2..9f1f57ddc5 100644 --- a/include/deal.II/lac/vector_space_vector.h +++ b/include/deal.II/lac/vector_space_vector.h @@ -77,7 +77,7 @@ namespace LinearAlgebra * communication pattern is used multiple times. This can be used to improve * performance. */ - virtual void Import(const ReadWriteVector &V, + virtual void import(const ReadWriteVector &V, VectorOperation::values operation, const CommunicationPatternBase *communication_pattern = NULL); diff --git a/source/lac/trilinos_epetra_vector.cc b/source/lac/trilinos_epetra_vector.cc index 6340d56fb9..b9ef42e3ae 100644 --- a/source/lac/trilinos_epetra_vector.cc +++ b/source/lac/trilinos_epetra_vector.cc @@ -33,14 +33,16 @@ namespace LinearAlgebra { Vector::Vector() : - vector(new Epetra_FEVector(Epetra_Map(0,1,0,Epetra_MpiComm(MPI_COMM_SELF)))) + vector(new Epetra_FEVector(Epetra_Map(0,1,0,Epetra_MpiComm(MPI_COMM_SELF)))), + epetra_comm_pattern(nullptr) {} Vector::Vector(const Vector &V) : - vector(new Epetra_FEVector(V.trilinos_vector())) + vector(new Epetra_FEVector(V.trilinos_vector())), + epetra_comm_pattern(nullptr) {} @@ -48,7 +50,8 @@ namespace LinearAlgebra Vector::Vector(const IndexSet ¶llel_partitioner, const MPI_Comm &communicator) : - vector(new Epetra_FEVector(parallel_partitioner.make_trilinos_map(communicator,false))) + vector(new Epetra_FEVector(parallel_partitioner.make_trilinos_map(communicator,false))), + epetra_comm_pattern(nullptr) {} @@ -96,91 +99,43 @@ namespace LinearAlgebra - Vector &Vector::operator= (const ::dealii::LinearAlgebra::ReadWriteVector &V) + void Vector::import(const ReadWriteVector &V, + VectorOperation::values operation, + const CommunicationPatternBase *communication_pattern) { - IndexSet elements_to_write(V.get_stored_elements()); - IndexSet locally_owned_elements(this->locally_owned_elements()); - IndexSet off_proc_elements(elements_to_write); - off_proc_elements.subtract_set(locally_owned_elements); - IndexSet on_proc_elements(elements_to_write); - on_proc_elements.subtract_set(off_proc_elements); - - // Write the elements that locally owned. - IndexSet::ElementIterator elem = on_proc_elements.begin(), - elem_end = on_proc_elements.end(); - for (; elem!=elem_end; ++elem) + // If no communication pattern is given, create one. Otherwsie, use the + // one given. + if (communication_pattern == nullptr) { - const TrilinosWrappers::types::int_type local_index = - vector->Map().LID(static_cast(*elem)); - (*vector)[0][local_index] = V[*elem]; - } - - - // Write the elements off-processor. Cannot use Import function of Trilinos - // because V is not a Trilinos object. The most straightforward to send - // the off-processor to the other processors is to AllGather. This way - // every processor has access to all the off-processor elements. The - // problem with this method is that if every element in V is - // off-processor, it is possible that the buffers used to receive the data - // has twice the global size of the current Vector. - // TODO This won't scale past a few 10,000's of processors. - - // Get recv_count and the size of the buffer needed to receive the data. - const Epetra_MpiComm *epetra_comm - = dynamic_cast(&(vector->Comm())); - MPI_Comm comm = epetra_comm->GetMpiComm(); - int send_buffer_size(off_proc_elements.n_elements()); - int n_procs(0); - MPI_Comm_size(comm, &n_procs); - std::vector recv_count(n_procs); - MPI_Allgather(&send_buffer_size, 1, MPI_INT, &recv_count, 1, - MPI_INT, comm); - std::vector displacement(n_procs,0); - for (int i=1; i send_index_buffer(send_buffer_size); - elem = off_proc_elements.begin(); - for (int i=0; i(vector->Comm()).Comm()); + } } - // Send the index of the off-processor elements - std::vector recv_index_buffer(recv_buffer_size); - MPI_Allgatherv(&send_index_buffer[0], send_buffer_size, DEAL_II_DOF_INDEX_MPI_TYPE, - &recv_index_buffer[0], &recv_count[0], &displacement[0], - DEAL_II_DOF_INDEX_MPI_TYPE, comm); - - // Write the values in the send buffer - std::vector send_value_buffer(send_buffer_size); - elem = off_proc_elements.begin(); - for (int i=0; i (communication_pattern)); + AssertThrow(epetra_comm_pattern != nullptr, + ExcMessage(std::string("The communication pattern is not of type ") + + "LinearAlgebra::EpetraWrappers::CommunicationPattern.")); } - // Send the value of the off-processor elements - std::vector recv_value_buffer(recv_buffer_size); - MPI_Allgatherv(&send_value_buffer[0], send_buffer_size, MPI_DOUBLE, - &recv_value_buffer[0], &recv_count[0], &displacement[0], MPI_DOUBLE, comm); - // Write the data in the vector - for (unsigned int i=0; iMap().LID(static_cast( - recv_index_buffer[i])); - (*vector)[0][local_index] = recv_value_buffer[i]; - } - } + Epetra_Import import(epetra_comm_pattern->get_epetra_import()); - return *this; + // The TargetMap and the SourceMap have their roles inverted. + Epetra_FEVector source_vector(import.TargetMap()); + double *values = source_vector.Values(); + std::copy(V.begin(), V.end(), values); + + if (operation==VectorOperation::insert) + vector->Export(source_vector, import, Insert); + else + vector->Export(source_vector, import, Add); } @@ -500,10 +455,10 @@ namespace LinearAlgebra out.setf(std::ios::fixed, std::ios::floatfield); if (across) - for (size_type i=0; iMyLength(); ++i) + for (int i=0; iMyLength(); ++i) out << val[i] << ' '; else - for (size_type i=0; iMyLength(); ++i) + for (int i=0; iMyLength(); ++i) out << val[i] << std::endl; out << std::endl; @@ -520,6 +475,16 @@ namespace LinearAlgebra + vector->MyLength()*(sizeof(double)+ sizeof(TrilinosWrappers::types::int_type)); } + + + + void Vector::create_epetra_comm_pattern(const IndexSet &source_index_set, + const MPI_Comm &mpi_comm) + { + source_stored_elements = source_index_set; + epetra_comm_pattern.reset(new CommunicationPattern(locally_owned_elements(), + source_index_set, mpi_comm)); + } } } -- 2.39.5