]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Make some functions const and use covariant return types. Add ComputationPattern
authorBruno Turcksin <bruno.turcksin@gmail.com>
Tue, 24 Nov 2015 22:45:10 +0000 (16:45 -0600)
committerBruno Turcksin <bruno.turcksin@gmail.com>
Wed, 23 Mar 2016 13:33:42 +0000 (09:33 -0400)
class.

13 files changed:
include/deal.II/lac/communication_pattern_base.h [new file with mode: 0644]
include/deal.II/lac/la_vector.h
include/deal.II/lac/la_vector.templates.h
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/trilinos_epetra_communication_pattern.h [new file with mode: 0644]
include/deal.II/lac/trilinos_epetra_vector.h
include/deal.II/lac/vector_space_vector.h
source/lac/CMakeLists.txt
source/lac/trilinos_epetra_communication_pattern.cc [new file with mode: 0644]
source/lac/trilinos_epetra_vector.cc
tests/trilinos/epetra_vector_01.cc
tests/trilinos/epetra_vector_02.cc

diff --git a/include/deal.II/lac/communication_pattern_base.h b/include/deal.II/lac/communication_pattern_base.h
new file mode 100644 (file)
index 0000000..9ef2852
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/mpi.h>
+
+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
index 09c70394af510d789c6ebfd7ec2260a7fafd32de..2ce6a28c82584c5f612cd45ae734efce49bcc1a0 100644 (file)
@@ -110,81 +110,81 @@ namespace LinearAlgebra
     /**
      * Multiply the entire vector by a fixed factor.
      */
-    virtual VectorSpaceVector<Number> &operator*= (const Number factor);
+    virtual Vector<Number> &operator*= (const Number factor) override;
 
     /**
      * Divide the entire vector by a fixed factor.
      */
-    virtual VectorSpaceVector<Number> &operator/= (const Number factor);
+    virtual Vector<Number> &operator/= (const Number factor) override;
 
     /**
      * Add the vector @p V to the present one.
      */
-    virtual VectorSpaceVector<Number> &operator+= (const VectorSpaceVector<Number> &V);
+    virtual Vector<Number> &operator+= (const VectorSpaceVector<Number> &V) override;
 
     /**
      * Substract the vector @p V from the present one.
      */
-    virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &V);
+    virtual Vector<Number> &operator-= (const VectorSpaceVector<Number> &V) override;
 
     /**
      * Return the scalar product of two vectors.
      */
-    virtual Number operator* (const VectorSpaceVector<Number> &V);
+    virtual Number operator* (const VectorSpaceVector<Number> &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. <tt>*this += a*V</tt>.
      */
-    virtual void add(const Number a, const VectorSpaceVector<Number> &V);
+    virtual void add(const Number a, const VectorSpaceVector<Number> &V) override;
 
     /**
      * Multiple addition of a multiple of a vector, i.e. <tt>*this +=
      * a*V+b*W</tt>.
      */
     virtual void add(const Number a, const VectorSpaceVector<Number> &V,
-                     const Number b, const VectorSpaceVector<Number> &W);
+                     const Number b, const VectorSpaceVector<Number> &W) override;
 
     /**
      * Scaling and simple addition of a multiple of a vector, i.e. <tt>*this =
      * s*(*this)+a*V</tt>.
      */
     virtual void sadd(const Number s, const Number a,
-                      const VectorSpaceVector<Number> &V);
+                      const VectorSpaceVector<Number> &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<Number> &scaling_factors);
+    virtual void scale(const VectorSpaceVector<Number> &scaling_factors) override;
 
     /**
      * Assignement <tt>*this = a*V</tt>.
      */
-    virtual void equ(const Number a, const VectorSpaceVector<Number> &V);
+    virtual void equ(const Number a, const VectorSpaceVector<Number> &V) override;
 
     /**
      * Return the l<sub>1</sub> norm of the vector (i.e., the sum of the
      * absolute values of all entries).
      */
-    virtual typename VectorSpaceVector<Number>::real_type l1_norm();
+    virtual typename VectorSpaceVector<Number>::real_type l1_norm() const override;
 
     /**
      * Return the l<sub>2</sub> norm of the vector (i.e., the square root of
      * the sum of the square of all entries among all processors).
      */
-    virtual typename VectorSpaceVector<Number>::real_type l2_norm();
+    virtual typename VectorSpaceVector<Number>::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<Number>::real_type linfty_norm();
+    virtual typename VectorSpaceVector<Number>::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<Number> &V,
-                               const VectorSpaceVector<Number> &W);
+                               const VectorSpaceVector<Number> &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<Number>::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<Number>::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 <typename Number>
   inline
-  const dealii::IndexSet Vector<Number>::locally_owned_elements() const
+  dealii::IndexSet Vector<Number>::locally_owned_elements() const
   {
-    return ReadWriteVector<Number>::get_stored_elements();
+    return IndexSet(ReadWriteVector<Number>::get_stored_elements());
   }
 
 
index 95e8465435667e747afae18e6348d37aea049db9..bc014e9daeb23e75d25b2cf904fe233755615d9a 100644 (file)
@@ -25,7 +25,7 @@ DEAL_II_NAMESPACE_OPEN
 namespace LinearAlgebra
 {
   template <typename Number>
-  VectorSpaceVector<Number> &Vector<Number>::operator*= (const Number factor)
+  Vector<Number> &Vector<Number>::operator*= (const Number factor)
   {
     AssertIsFinite(factor);
     for (unsigned int i=0; i<this->size(); ++i)
@@ -37,7 +37,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  VectorSpaceVector<Number> &Vector<Number>::operator/= (const Number factor)
+  Vector<Number> &Vector<Number>::operator/= (const Number factor)
   {
     AssertIsFinite(factor);
     Assert(factor!=Number(0.), ExcZero());
@@ -49,7 +49,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  VectorSpaceVector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
+  Vector<Number> &Vector<Number>::operator+= (const VectorSpaceVector<Number> &V)
   {
     // Check that casting will work.
     Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
@@ -70,7 +70,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  VectorSpaceVector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
+  Vector<Number> &Vector<Number>::operator-= (const VectorSpaceVector<Number> &V)
   {
     // Check that casting will work.
     Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL, ExcVectorTypeNotCompatible());
@@ -91,7 +91,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V)
+  Number Vector<Number>::operator* (const VectorSpaceVector<Number> &V) const
   {
     // Check that casting will work.
     Assert(dynamic_cast<const Vector<Number>*>(&V)!=NULL,
@@ -221,7 +221,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm()
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm() const
   {
     Assert (this->size(), ExcEmptyObject());
 
@@ -235,7 +235,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm()
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm() const
   {
     Assert (this->size(), ExcEmptyObject());
 
@@ -277,7 +277,7 @@ namespace LinearAlgebra
 
 
   template <typename Number>
-  typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm()
+  typename VectorSpaceVector<Number>::real_type Vector<Number>::linfty_norm() const
   {
     typename ReadWriteVector<Number>::real_type norm = 0.;
     for (unsigned int i=0; i<this->size(); ++i)
@@ -303,7 +303,7 @@ namespace LinearAlgebra
   template <typename Number>
   typename VectorSpaceVector<Number>::real_type Vector<Number>::l1_norm_recursive(
     unsigned int i,
-    unsigned int j)
+    unsigned int j) const
   {
     Assert(j>=i, ExcInternalError());
     typename ReadWriteVector<Number>::real_type norm = 0.;
@@ -324,7 +324,7 @@ namespace LinearAlgebra
   template <typename Number>
   typename VectorSpaceVector<Number>::real_type Vector<Number>::l2_norm_squared_recursive(
     unsigned int i,
-    unsigned int j)
+    unsigned int j) const
   {
     Assert(j>=i, ExcInternalError());
 
index a2574caa0729dd4ebfbca4c4af6c2a7ab0b7021b..f39e50c8de70fb673d37b9a496f9e19f545034a9 100644 (file)
 
 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<Number> &
-    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<Number> &
-    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<Number> &
-    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<Number> &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<IndexSet> source_stored_elements;
+
+    /**
+     *
+     */
+    std::shared_ptr<CommunicationPatternBase> comm_pattern;
+
     /**
      * Pointer to the array of local elements of this vector.
      */
@@ -471,7 +511,9 @@ namespace LinearAlgebra
   inline
   ReadWriteVector<Number>::ReadWriteVector ()
     :
-    val (NULL)
+    source_stored_elements(nullptr),
+    comm_pattern(nullptr),
+    val(nullptr)
   {}
 
 
@@ -481,7 +523,9 @@ namespace LinearAlgebra
   ReadWriteVector<Number>::ReadWriteVector (const ReadWriteVector<Number> &v)
     :
     Subscriptor(),
-    val (NULL)
+    source_stored_elements(nullptr),
+    comm_pattern(nullptr),
+    val(nullptr)
   {
     this->operator=(v);
   }
@@ -492,7 +536,9 @@ namespace LinearAlgebra
   inline
   ReadWriteVector<Number>::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<Number>::ReadWriteVector (const IndexSet &locally_stored_indices)
     :
-    val (NULL)
+    source_stored_elements(nullptr),
+    comm_pattern(nullptr),
+    val(nullptr)
   {
     reinit (locally_stored_indices);
   }
index 936574e918de2815171155ae42c43e272b83e4fb..0e42ecf04bcdb220212bbdb684ecf4e3272b0b7d 100644 (file)
@@ -1,6 +1,6 @@
 // ---------------------------------------------------------------------
 //
-// Copyright (C) 20152016 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 <deal.II/lac/trilinos_epetra_vector.h>
 
 #ifdef DEAL_II_WITH_TRILINOS
-#  include <Epetra_Import.h>
+#  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
+#  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 <typename Number>
-  ReadWriteVector<Number> &
-  ReadWriteVector<Number>::operator = (const PETScWrappers::MPI::Vector &petsc_vec)
+  void
+  ReadWriteVector<Number>::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<const Vec &>(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 <typename Number>
-  ReadWriteVector<Number> &
-  ReadWriteVector<Number>::import(const Epetra_MultiVector &multivector,
-                                  const IndexSet           &locally_owned_elements)
+  void
+  ReadWriteVector<Number>::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<size_type> trilinos_indices, readwrite_indices;
-    locally_owned_elements.fill_index_vector(trilinos_indices);
-    stored_elements.fill_index_vector(readwrite_indices);
-    std::vector<size_type> intersection(trilinos_indices.size());
-    std::vector<size_type>::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<size_t> trilinos_position(intersection_size);
-    unsigned int j=0;
-    for (size_t i=0; i<intersection_size; ++i)
+    const EpetraWrappers::CommunicationPattern *epetra_comm_pattern = nullptr;
+
+    if (communication_pattern == nullptr)
       {
-        while (intersection[i]!=trilinos_indices[j])
-          ++j;
-        trilinos_position[i] = j;
+        if (source_elements == *source_stored_elements)
+          {
+            epetra_comm_pattern =
+              dynamic_cast<const EpetraWrappers::CommunicationPattern *> (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<const EpetraWrappers::CommunicationPattern *> (comm_pattern.get());
+      }
+    else
+      {
+        epetra_comm_pattern =
+          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."));
       }
-    double *values = multivector.Values();
-    for (size_t i=0; i<trilinos_position.size(); ++i)
-      val[global_to_local(intersection[i])] = values[trilinos_position[i]];
 
+    Epetra_Import import(epetra_comm_pattern->get_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<size_type> 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.size(); ++i)
-          off_proc_array[i] = static_cast<TrilinosWrappers::types::int_type> (
-                                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<off_proc_indices.size(); ++i)
-          val[global_to_local(off_proc_indices[i])] = values[i];
+        double *values = target_vector.Values();
+        for (int i=0; i<target_vector.MyLength(); ++i)
+          val[i] = values[i];
+      }
+    else
+      {
+        target_vector.Import(multivector, import, Add);
+        double *values = target_vector.Values();
+        for (int i=0; i<target_vector.MyLength(); ++i)
+          val[i] += values[i];
       }
-#endif
-
-    // return a pointer to this object per normal c++ operator overloading
-    // semantics
-    return *this;
   }
-#endif
 
 
   template <typename Number>
-  ReadWriteVector<Number> &
-  ReadWriteVector<Number>::operator = (const TrilinosWrappers::MPI::Vector &trilinos_vec)
+  void
+  ReadWriteVector<Number>::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 <typename Number>
-  ReadWriteVector<Number> &
-  ReadWriteVector<Number>::operator = (const LinearAlgebra::EpetraWrappers::Vector &trilinos_vec)
+  void
+  ReadWriteVector<Number>::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 <typename Number>
@@ -296,6 +307,18 @@ namespace LinearAlgebra
     out.flags (old_flags);
     out.precision(old_precision);
   }
+
+
+
+  template <typename Number>
+  void
+  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));
+  }
 } // 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 (file)
index 0000000..d25f918
--- /dev/null
@@ -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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/std_cxx11/shared_ptr.h>
+#include <deal.II/lac/communication_pattern_base.h>
+#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<const MPI_Comm> comm;
+
+      /**
+       * Shared pointer to the Epetra_Import object used.
+       */
+      std_cxx11::shared_ptr<Epetra_Import> import;
+    };
+  } // end of namespace EpetraWrappers
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
+
+#endif
index 61301b6f6b9a1278c868bb340acd47c251766a29..527a9af2548cc1b9f3b3aa5aa7b97abdd941970e 100644 (file)
@@ -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 &parallel_partitioner,
-             const MPI_Comm &communicator = MPI_COMM_WORLD);
+      explicit Vector(const IndexSet &parallel_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 &parallel_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<double> &operator*= (const double factor);
+      virtual Vector &operator*= (const double factor);
 
       /**
        * Divide the entire vector by a fixed factor.
        */
-      virtual VectorSpaceVector<double> &operator/= (const double factor);
+      virtual Vector &operator/= (const double factor);
 
       /**
        * Add the vector @p V to the present one.
        */
-      virtual VectorSpaceVector<double> &operator+= (const VectorSpaceVector<double> &V);
+      virtual Vector &operator+= (const VectorSpaceVector<double> &V);
 
       /**
        * Substract the vector @p V from the present one.
        */
-      virtual VectorSpaceVector<double> &operator-= (const VectorSpaceVector<double> &V);
+      virtual Vector &operator-= (const VectorSpaceVector<double> &V);
 
       /**
        * Return the scalar product of two vectors. The vectors need to have the
        * same layout.
        */
-      virtual double operator* (const VectorSpaceVector<double> &V);
+      virtual double operator* (const VectorSpaceVector<double> &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 l<sub>1</sub> 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 l<sub>2</sub> 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
index ec008f89d3271c4f9d065546b1fddc84dbae6b2c..26c0b584e2fe2d147eaa74d6f5ea4ba6ef8bea5d 100644 (file)
@@ -69,10 +69,22 @@ namespace LinearAlgebra
      */
     virtual VectorSpaceVector<Number> &operator-= (const VectorSpaceVector<Number> &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<Number> &V,
+                        VectorOperation::values operation,
+                        const CommunicationPatternBase *communication_pattern = NULL);
+
     /**
      * Return the scalar product of two vectors.
      */
-    virtual Number operator* (const VectorSpaceVector<Number> &V) = 0;
+    virtual Number operator* (const VectorSpaceVector<Number> &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 l<sub>1</sub> 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 l<sub>2</sub> 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.
index b131f2746f8021f8a3548f3b0c8a55fc14f5ffed..9dc145c2c5943827853d36622409fa4998343262 100644 (file)
@@ -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 (file)
index 0000000..ef93fc9
--- /dev/null
@@ -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 <deal.II/lac/trilinos_epetra_communication_pattern.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#ifdef DEAL_II_WITH_MPI
+
+#include <deal.II/base/index_set.h>
+#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<const MPI_Comm>(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
index 4bf7b32ea90c90b62093bc02520fd45db66e2a83..6340d56fb989cd1c4822eb6f70ee323c335418e1 100644 (file)
@@ -16,6 +16,7 @@
 #include <deal.II/lac/trilinos_epetra_vector.h>
 
 #ifdef DEAL_II_WITH_TRILINOS
+
 #ifdef DEAL_II_WITH_MPI
 
 #include <deal.II/base/index_set.h>
@@ -184,7 +185,7 @@ namespace LinearAlgebra
 
 
 
-    VectorSpaceVector<double> &Vector::operator*= (const double factor)
+    Vector &Vector::operator*= (const double factor)
     {
       AssertIsFinite(factor);
       vector->Scale(factor);
@@ -194,7 +195,7 @@ namespace LinearAlgebra
 
 
 
-    VectorSpaceVector<double> &Vector::operator/= (const double factor)
+    Vector &Vector::operator/= (const double factor)
     {
       AssertIsFinite(factor);
       Assert(factor!=0., ExcZero());
@@ -205,7 +206,7 @@ namespace LinearAlgebra
 
 
 
-    VectorSpaceVector<double> &Vector::operator+= (const VectorSpaceVector<double> &V)
+    Vector &Vector::operator+= (const VectorSpaceVector<double> &V)
     {
       // Check that casting will work.
       Assert(dynamic_cast<const Vector *>(&V)!=NULL, ExcVectorTypeNotCompatible());
@@ -247,7 +248,7 @@ namespace LinearAlgebra
 
 
 
-    VectorSpaceVector<double> &Vector::operator-= (const VectorSpaceVector<double> &V)
+    Vector &Vector::operator-= (const VectorSpaceVector<double> &V)
     {
       this->add(-1.,V);
 
@@ -256,7 +257,7 @@ namespace LinearAlgebra
 
 
 
-    double Vector::operator* (const VectorSpaceVector<double> &V)
+    double Vector::operator* (const VectorSpaceVector<double> &V) const
     {
       // Check that casting will work.
       Assert(dynamic_cast<const Vector *>(&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<const Epetra_Map *>(&(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; i<size(); ++i)
+        for (size_type i=0; i<vector->MyLength(); ++i)
           out << val[i] << ' ';
       else
-        for (size_type i=0; i<size(); ++i)
+        for (size_type i=0; i<vector->MyLength(); ++i)
           out << val[i] << std::endl;
       out << std::endl;
 
index c8fe4ee4226ade1278c5d28a3e7b458afcf162fc..3059910224c34fb9afb480e443aa3e63af8325c6 100644 (file)
@@ -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(),
index 746de05288a76b868c386af84490d0680be38063..687038a39397dc98fc3ac61545f0a8702b165049 100644 (file)
@@ -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)

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.