]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add import function to LinearAlgebra::TrilinosWrappers::Vector.
authorBruno Turcksin <bruno.turcksin@gmail.com>
Fri, 12 Feb 2016 19:13:36 +0000 (14:13 -0500)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 23 Mar 2016 13:33:42 +0000 (09:33 -0400)
include/deal.II/lac/communication_pattern_base.h
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/trilinos_epetra_vector.h
include/deal.II/lac/vector_space_vector.h
source/lac/trilinos_epetra_vector.cc

index 9ef2852291eaf17192bf9dcdcf9168abf2306fe4..e63a2d25fd111de64daafe6bb7fe354e7df8b69c 100644 (file)
@@ -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.
    */
index f39e50c8de70fb673d37b9a496f9e19f545034a9..cb61c2c9a705a3055453a2c251955497a9210394 100644 (file)
@@ -30,6 +30,7 @@
 #include <iomanip>
 
 #ifdef DEAL_II_WITH_TRILINOS
+#include <deal.II/lac/trilinos_epetra_communication_pattern.h>
 #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<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;
 
@@ -511,7 +517,6 @@ namespace LinearAlgebra
   inline
   ReadWriteVector<Number>::ReadWriteVector ()
     :
-    source_stored_elements(nullptr),
     comm_pattern(nullptr),
     val(nullptr)
   {}
@@ -523,7 +528,6 @@ namespace LinearAlgebra
   ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
     :
     Subscriptor(),
-    source_stored_elements(nullptr),
     comm_pattern(nullptr),
     val(nullptr)
   {
@@ -536,7 +540,6 @@ namespace LinearAlgebra
   inline
   ReadWriteVector<Number>::ReadWriteVector (const size_type size)
     :
-    source_stored_elements(nullptr),
     comm_pattern(nullptr),
     val(nullptr)
   {
@@ -549,7 +552,6 @@ namespace LinearAlgebra
   inline
   ReadWriteVector<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
     :
-    source_stored_elements(nullptr),
     comm_pattern(nullptr),
     val(nullptr)
   {
index 0e42ecf04bcdb220212bbdb684ecf4e3272b0b7d..e3a7ea8a56d7d30c8a62a2f991c996aba053e69a 100644 (file)
@@ -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<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."));
@@ -310,15 +314,21 @@ namespace LinearAlgebra
 
 
 
+#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
 
 
index 527a9af2548cc1b9f3b3aa5aa7b97abdd941970e..67fd2efe714bf08ec6302e68735daa7743479d9f 100644 (file)
@@ -27,6 +27,7 @@
 #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"
 
@@ -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<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.
@@ -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<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;
     };
   }
 }
index 26c0b584e2fe2d147eaa74d6f5ea4ba6ef8bea5d..9f1f57ddc509ad0ba3c36f8612ee0b70d5ac5a93 100644 (file)
@@ -77,7 +77,7 @@ namespace LinearAlgebra
      * 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);
 
index 6340d56fb989cd1c4822eb6f70ee323c335418e1..b9ef42e3aeb1f056487424f44084a2e3b0ee9c6f 100644 (file)
@@ -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 &parallel_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<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);
     }
 
 
@@ -500,10 +455,10 @@ namespace LinearAlgebra
         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;
 
@@ -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));
+    }
   }
 }
 

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.