/**
* 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.
*/
#include <iomanip>
#ifdef DEAL_II_WITH_TRILINOS
+#include <deal.II/lac/trilinos_epetra_communication_pattern.h>
#include "Epetra_MultiVector.h"
#endif
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
*/
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.
IndexSet stored_elements;
/**
- * Indices of the elements stored on a VectorSpaceVector
+ * IndexSet of the elements of the last imported vector;
*/
- std::shared_ptr<IndexSet> source_stored_elements;
+ IndexSet source_stored_elements;
/**
- *
+ * CommunicationPattern for the communication between the
+ * source_stored_elements IndexSet and the current vector.
*/
std::shared_ptr<CommunicationPatternBase> comm_pattern;
inline
ReadWriteVector<Number>::ReadWriteVector ()
:
- source_stored_elements(nullptr),
comm_pattern(nullptr),
val(nullptr)
{}
ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
:
Subscriptor(),
- source_stored_elements(nullptr),
comm_pattern(nullptr),
val(nullptr)
{
inline
ReadWriteVector<Number>::ReadWriteVector (const size_type size)
:
- source_stored_elements(nullptr),
comm_pattern(nullptr),
val(nullptr)
{
inline
ReadWriteVector<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
:
- source_stored_elements(nullptr),
comm_pattern(nullptr),
val(nullptr)
{
this->operator = (Number());
// reset the communication patter
- source_stored_elements = nullptr;
+ source_stored_elements.clear();
comm_pattern = nullptr;
}
this->operator= (Number());
// reset the communication patter
- source_stored_elements = nullptr;
+ source_stored_elements.clear();
comm_pattern = nullptr;
}
this->operator= (Number());
// reset the communication patter
- source_stored_elements = nullptr;
+ source_stored_elements.clear();
comm_pattern = nullptr;
}
const MPI_Comm &mpi_comm,
const CommunicationPatternBase *communication_pattern)
{
- const EpetraWrappers::CommunicationPattern *epetra_comm_pattern = nullptr;
+ std_cxx11::shared_ptr<const EpetraWrappers::CommunicationPattern> 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<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get());
+ epetra_comm_pattern.reset(
+ dynamic_cast<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get()));
if (epetra_comm_pattern == nullptr)
- create_epetra_comm_pattern(source_elements, mpi_comm);
+ epetra_comm_pattern = std::make_shared<const EpetraWrappers::CommunicationPattern>(
+ create_epetra_comm_pattern(source_elements, mpi_comm));
}
else
- create_epetra_comm_pattern(source_elements, mpi_comm);
-
- epetra_comm_pattern =
- dynamic_cast<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get());
+ epetra_comm_pattern = std::make_shared<const EpetraWrappers::CommunicationPattern>(
+ create_epetra_comm_pattern(source_elements, mpi_comm));
}
else
{
- epetra_comm_pattern =
- dynamic_cast<const EpetraWrappers::CommunicationPattern *> (communication_pattern);
+ epetra_comm_pattern.reset(
+ dynamic_cast<const EpetraWrappers::CommunicationPattern *> (communication_pattern));
AssertThrow(epetra_comm_pattern != nullptr,
ExcMessage(std::string("The communication pattern is not of type ") +
"LinearAlgebra::EpetraWrappers::CommunicationPattern."));
+#ifdef DEAL_II_WITH_TRILINOS
template <typename Number>
- void
+ EpetraWrappers::CommunicationPattern
ReadWriteVector<Number>::create_epetra_comm_pattern(const IndexSet &source_index_set,
const MPI_Comm &mpi_comm)
{
- source_stored_elements = std::make_shared<IndexSet>(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<EpetraWrappers::CommunicationPattern>(
+ epetra_comm_pattern);
+
+ return epetra_comm_pattern;
}
+#endif
} // end of namespace LinearAlgebra
#include <deal.II/base/std_cxx11/shared_ptr.h>
#include <deal.II/lac/read_write_vector.h>
#include <deal.II/lac/vector_space_vector.h>
+#include <memory>
#include "mpi.h"
#include "Epetra_FEVector.h"
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<double> &V);
+ virtual void import(const ReadWriteVector<double> &V,
+ VectorOperation::values operation,
+ const CommunicationPatternBase *communication_pattern = NULL)
+ override;
/**
* Multiply the entire vector by a fixed factor.
<< " 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<Epetra_FEVector> 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<const CommunicationPattern> epetra_comm_pattern;
};
}
}
* communication pattern is used multiple times. This can be used to improve
* performance.
*/
- virtual void Import(const ReadWriteVector<Number> &V,
+ virtual void import(const ReadWriteVector<Number> &V,
VectorOperation::values operation,
const CommunicationPatternBase *communication_pattern = NULL);
{
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)
{}
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)
{}
- Vector &Vector::operator= (const ::dealii::LinearAlgebra::ReadWriteVector<double> &V)
+ void Vector::import(const ReadWriteVector<double> &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<TrilinosWrappers::types::int_type>(*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<const Epetra_MpiComm *>(&(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<int> recv_count(n_procs);
- MPI_Allgather(&send_buffer_size, 1, MPI_INT, &recv_count, 1,
- MPI_INT, comm);
- std::vector<int> displacement(n_procs,0);
- for (int i=1; i<n_procs; ++i)
- displacement[i] = displacement[i-1] + recv_count[i-1];
- const unsigned int recv_buffer_size = displacement[n_procs-1] +
- recv_count[n_procs-1];
-
- // Write the indices in the send buffer
- std::vector<types::global_dof_index> send_index_buffer(send_buffer_size);
- elem = off_proc_elements.begin();
- for (int i=0; i<send_buffer_size; ++i)
- {
- send_index_buffer[i] = *elem;
- ++elem;
+ // The first time import is called, a communication pattern is created.
+ // Check if the communication pattern already exists and if it can be
+ // reused.
+ if (source_stored_elements == V.get_stored_elements())
+ {
+ create_epetra_comm_pattern(V.get_stored_elements(),
+ dynamic_cast<const Epetra_MpiComm &>(vector->Comm()).Comm());
+ }
}
- // Send the index of the off-processor elements
- std::vector<types::global_dof_index> 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<double> send_value_buffer(send_buffer_size);
- elem = off_proc_elements.begin();
- for (int i=0; i<send_buffer_size; ++i)
+ else
{
- send_value_buffer[i] = V[*elem];
- ++elem;
+ epetra_comm_pattern.reset(
+ dynamic_cast<const CommunicationPattern *> (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<double> 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; i<recv_buffer_size; ++i)
- {
- if (locally_owned_elements.is_element(recv_index_buffer[i]))
- {
- const TrilinosWrappers::types::int_type local_index =
- vector->Map().LID(static_cast<TrilinosWrappers::types::int_type>(
- 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);
}
out.setf(std::ios::fixed, std::ios::floatfield);
if (across)
- for (size_type i=0; i<vector->MyLength(); ++i)
+ for (int i=0; i<vector->MyLength(); ++i)
out << val[i] << ' ';
else
- for (size_type i=0; i<vector->MyLength(); ++i)
+ for (int i=0; i<vector->MyLength(); ++i)
out << val[i] << std::endl;
out << std::endl;
+ 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));
+ }
}
}
projects / dealii.git / commitdiff