]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement TpetraWrappers::Vector
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Mon, 27 Aug 2018 15:14:38 +0000 (17:14 +0200)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Wed, 20 Feb 2019 12:46:19 +0000 (13:46 +0100)
39 files changed:
cmake/config/template-arguments.in
cmake/configure/configure_2_trilinos.cmake
doc/external-libs/trilinos.html
include/deal.II/base/index_set.h
include/deal.II/base/utilities.h
include/deal.II/dofs/dof_accessor.templates.h
include/deal.II/fe/fe_tools_extrapolate.templates.h
include/deal.II/fe/fe_tools_interpolate.templates.h
include/deal.II/lac/read_write_vector.h
include/deal.II/lac/read_write_vector.templates.h
include/deal.II/lac/trilinos_sparse_matrix.h
include/deal.II/lac/trilinos_tpetra_communication_pattern.h [new file with mode: 0644]
include/deal.II/lac/trilinos_tpetra_vector.h [new file with mode: 0644]
include/deal.II/lac/vector_element_access.h
include/deal.II/multigrid/mg_transfer.h
include/deal.II/numerics/vector_tools.templates.h
source/base/index_set.cc
source/base/time_stepping.cc
source/base/utilities.cc
source/dofs/dof_accessor_get.cc
source/dofs/dof_accessor_set.cc
source/fe/fe_values.cc
source/fe/mapping_fe_field.cc
source/lac/CMakeLists.txt
source/lac/solver.cc
source/lac/trilinos_sparse_matrix.cc
source/lac/trilinos_tpetra_communication_pattern.cc [new file with mode: 0644]
source/lac/trilinos_tpetra_vector.cc [new file with mode: 0644]
source/lac/vector_memory.cc
source/multigrid/mg_base.cc
source/numerics/data_out_dof_data_codim.cc
source/numerics/derivative_approximation.cc
tests/trilinos/tpetra_vector_01.cc [new file with mode: 0644]
tests/trilinos/tpetra_vector_01.mpirun=2.output [new file with mode: 0644]
tests/trilinos/tpetra_vector_02.cc [new file with mode: 0644]
tests/trilinos/tpetra_vector_02.mpirun=2.output [new file with mode: 0644]
tests/trilinos/tpetra_vector_03.cc [new file with mode: 0644]
tests/trilinos/tpetra_vector_03.mpirun=2.output [new file with mode: 0644]
tests/trilinos/tpetra_vector_03.mpirun=4.output [new file with mode: 0644]

index 551fdd08d3cc8ca8b1fd0ea962526bd74997eea0..17a9f2ba94776fe00e3ce9271db04d6e265c25b5 100644 (file)
@@ -95,6 +95,7 @@ VECTOR_TYPES := {   Vector<double>;
 
                     @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                     @DEAL_II_EXPAND_EPETRA_VECTOR@;
+                    @DEAL_II_EXPAND_TPETRA_VECTOR@;
                     @DEAL_II_EXPAND_PETSC_MPI_VECTOR@;
 
                     @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@;
@@ -120,6 +121,7 @@ REAL_VECTOR_TYPES  := {  Vector<double>;
 
                          @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                          @DEAL_II_EXPAND_EPETRA_VECTOR@;
+                         @DEAL_II_EXPAND_TPETRA_VECTOR@;
                          @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@;
 
                          @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@;
@@ -138,6 +140,7 @@ REAL_NONBLOCK_VECTORS := { Vector<double>;
 
                            @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                            @DEAL_II_EXPAND_EPETRA_VECTOR@;
+                           @DEAL_II_EXPAND_TPETRA_VECTOR@;
                            @DEAL_II_EXPAND_PETSC_MPI_VECTOR_REAL@;
                          }
 
@@ -145,6 +148,7 @@ REAL_NONBLOCK_VECTORS := { Vector<double>;
 EXTERNAL_PARALLEL_VECTORS := { @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                                @DEAL_II_EXPAND_TRILINOS_MPI_BLOCKVECTOR@;
                                @DEAL_II_EXPAND_EPETRA_VECTOR@;
+                               @DEAL_II_EXPAND_TPETRA_VECTOR@;
                                @DEAL_II_EXPAND_PETSC_MPI_VECTOR@;
                                @DEAL_II_EXPAND_PETSC_MPI_BLOCKVECTOR@
                              }
@@ -163,6 +167,7 @@ VECTORS_WITH_MATRIX := { Vector<double>;
 
                          @DEAL_II_EXPAND_TRILINOS_MPI_VECTOR@;
                          @DEAL_II_EXPAND_EPETRA_VECTOR@;
+                         @DEAL_II_EXPAND_TPETRA_VECTOR@;
                          @DEAL_II_EXPAND_PETSC_MPI_VECTOR@;
                        }
 
index c53f105c0995ba699b105cc9b56bee39fa71d037..abe66905d82c2f73b9fb790afa1c9c7624247a75 100644 (file)
@@ -42,7 +42,7 @@ MACRO(FEATURE_TRILINOS_FIND_EXTERNAL var)
       )
 
     FOREACH(_module
-        Amesos Epetra Ifpack AztecOO Teuchos ML MueLu
+        Amesos Epetra Ifpack AztecOO Teuchos Tpetra ML MueLu
       )
       ITEM_MATCHES(_module_found ${_module} ${Trilinos_PACKAGE_LIST})
       IF(_module_found)
@@ -227,6 +227,7 @@ MACRO(FEATURE_TRILINOS_CONFIGURE_EXTERNAL)
   SET(DEAL_II_EXPAND_TRILINOS_MPI_VECTOR "TrilinosWrappers::MPI::Vector")
   IF (TRILINOS_WITH_MPI)
     SET(DEAL_II_EXPAND_EPETRA_VECTOR "LinearAlgebra::EpetraWrappers::Vector")
+    SET(DEAL_II_EXPAND_TPETRA_VECTOR "LinearAlgebra::TpetraWrappers::Vector")
   ENDIF()
   IF(${DEAL_II_TRILINOS_WITH_SACADO})
     # Note: Only CMake 3.0 and greater support line continuation with the "\" character
index 9579c9a62a0a6c314599e2216c646733f7c4a48e..d950f55b11c096c65eb96c8209afb042295b066d 100644 (file)
@@ -67,6 +67,7 @@
         <li> ROL (optional),
         <li> Sacado (optional),
         <li> Teuchos,
+        <li> Tpetra,
         <li> Zoltan (optional).
       </ul>
 
@@ -92,6 +93,7 @@
     -DTrilinos_ENABLE_MueLu=ON                       \
     -DTrilinos_ENABLE_ML=ON                          \
     -DTrilinos_ENABLE_ROL=ON                         \
+    -DTrilinos_ENABLE_Tpetra=ON                      \
     -DTrilinos_ENABLE_Zoltan=ON                      \
     -DTrilinos_VERBOSE_CONFIGURE=OFF                 \
     -DTPL_ENABLE_MPI=ON                              \
index fec38b883f22fdb178578a1b1158a94032d747d3..2e21fab6449883be1a565c54663e8b31e7d4b5ae 100644 (file)
@@ -31,6 +31,7 @@
 
 #ifdef DEAL_II_WITH_TRILINOS
 #  include <Epetra_Map.h>
+#  include <Tpetra_Map.hpp>
 #endif
 
 #if defined(DEAL_II_WITH_MPI) || defined(DEAL_II_WITH_PETSC)
@@ -446,6 +447,10 @@ public:
   Epetra_Map
   make_trilinos_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
                     const bool      overlapping  = false) const;
+
+  Tpetra::Map<>
+  make_tpetra_map(const MPI_Comm &communicator = MPI_COMM_WORLD,
+                  const bool      overlapping  = false) const;
 #endif
 
 
index f49c3021321b07040b3bb5086cc86a13646c2203..71d05faf24ea5beb844768ab26664118bd38921d 100644 (file)
@@ -888,6 +888,9 @@ namespace Utilities
     const Epetra_Comm &
     comm_self();
 
+    const Teuchos::RCP<const Teuchos::Comm<int>> &
+    tpetra_comm_self();
+
     /**
      * Given a communicator, duplicate it. If the given communicator is
      * serial, that means to just return a copy of itself. On the other hand,
index a7c898c9537f134d5962a7a2c2b79a231721689d..73b77c07a57d79e625ca721bc594916578c600b9 100644 (file)
@@ -1534,6 +1534,33 @@ namespace internal
 
 
 
+      template <typename ForwardIterator>
+      static void
+      extract_subvector_to(const LinearAlgebra::TpetraWrappers::Vector &values,
+                           const types::global_dof_index *cache_begin,
+                           const types::global_dof_index *cache_end,
+                           ForwardIterator                local_values_begin)
+      {
+        std::vector<unsigned int> sorted_indices_pos =
+          sort_indices(cache_begin, cache_end);
+        const unsigned int cache_size = cache_end - cache_begin;
+        std::vector<types::global_dof_index> cache_indices(cache_size);
+        for (unsigned int i = 0; i < cache_size; ++i)
+          cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
+
+        IndexSet index_set(cache_indices.back() + 1);
+        index_set.add_indices(cache_indices.begin(), cache_indices.end());
+        index_set.compress();
+        LinearAlgebra::ReadWriteVector<double> read_write_vector(index_set);
+        read_write_vector.import(values, VectorOperation::insert);
+
+        // Copy the elements from read_write_vector and reorder them.
+        for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
+          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
+      }
+
+
+
       template <typename ForwardIterator>
       static void
       extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values,
index f2acc692d7aec3c6ba6a44ad80d984d5cc1f7fcc..850a4c6216cab919e741f89a0bfc6152d27403d6 100644 (file)
@@ -1581,6 +1581,21 @@ namespace FETools
 
 
 #  ifdef DEAL_II_WITH_MPI
+    template <int dim, int spacedim>
+    void
+    reinit_distributed(const DoFHandler<dim, spacedim> &      dh,
+                       LinearAlgebra::TpetraWrappers::Vector &vector)
+    {
+      const parallel::distributed::Triangulation<dim, spacedim> *parallel_tria =
+        dynamic_cast<
+          const parallel::distributed::Triangulation<dim, spacedim> *>(
+          &dh.get_triangulation());
+      Assert(parallel_tria != nullptr, ExcNotImplemented());
+
+      const IndexSet &locally_owned_dofs = dh.locally_owned_dofs();
+      vector.reinit(locally_owned_dofs, parallel_tria->get_communicator());
+    }
+
     template <int dim, int spacedim>
     void
     reinit_distributed(const DoFHandler<dim, spacedim> &      dh,
index 82460f14bb248c497c3bff1a6622b0713c5d4b25..e2e9b417fe02a87316d1935c459157ef99c52254 100644 (file)
@@ -46,6 +46,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
@@ -511,6 +512,21 @@ namespace FETools
     {
       AssertThrow(false, ExcNotImplemented());
     }
+
+    template <int dim, int spacedim>
+    void
+    back_interpolate(
+      const DoFHandler<dim, spacedim> &,
+      const AffineConstraints<
+        typename LinearAlgebra::TpetraWrappers::Vector::value_type> &,
+      const LinearAlgebra::TpetraWrappers::Vector &,
+      const DoFHandler<dim, spacedim> &,
+      const AffineConstraints<
+        typename LinearAlgebra::TpetraWrappers::Vector::value_type> &,
+      LinearAlgebra::TpetraWrappers::Vector &)
+    {
+      AssertThrow(false, ExcNotImplemented());
+    }
 #endif
 
 
index 9e67f981787d07044c2b59161d1b9bad600bd79a..70d89ade7a1c88a3d3a2fb0de2c762f9fc07a567 100644 (file)
@@ -35,6 +35,7 @@
 #ifdef DEAL_II_WITH_TRILINOS
 #  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
 #  include <deal.II/lac/trilinos_epetra_vector.h>
+#  include <deal.II/lac/trilinos_tpetra_vector.h>
 
 #  include <Epetra_MultiVector.h>
 
@@ -334,6 +335,21 @@ namespace LinearAlgebra
                std::shared_ptr<const CommunicationPatternBase>());
 
 #  ifdef DEAL_II_WITH_MPI
+    /**
+     * Imports all the elements present in the vector's IndexSet from the input
+     * vector @p tpetra_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 TpetraWrappers::Vector &tpetra_vec,
+           VectorOperation::values       operation,
+           const std::shared_ptr<const CommunicationPatternBase>
+             &communication_pattern =
+               std::shared_ptr<const CommunicationPatternBase>());
+
     /**
      * Imports all the elements present in the vector's IndexSet from the input
      * vector @p epetra_vec. VectorOperation::values @p operation is used to
@@ -593,6 +609,19 @@ namespace LinearAlgebra
 
   protected:
 #ifdef DEAL_II_WITH_TRILINOS
+    /**
+     * Import all the elements present in the vector's IndexSet from the input
+     * vector @p tpetra_vector. This is an helper function and it should not be
+     * used directly.
+     */
+    void
+    import(const Tpetra::Vector<> &tpetra_vector,
+           const IndexSet &        locally_owned_elements,
+           VectorOperation::values operation,
+           const MPI_Comm &        mpi_comm,
+           const std::shared_ptr<const CommunicationPatternBase>
+             &communication_pattern);
+
     /**
      * Import all the elements present in the vector's IndexSet from the input
      * vector @p multivector. This is an helper function and it should not be
@@ -627,7 +656,15 @@ namespace LinearAlgebra
 
 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
     /**
-     * Return a EpetraWrappers::Communication pattern and store it for future
+     * Return a TpetraWrappers::CommunicationPattern and store it for future
+     * use.
+     */
+    TpetraWrappers::CommunicationPattern
+    create_tpetra_comm_pattern(const IndexSet &source_index_set,
+                               const MPI_Comm &mpi_comm);
+
+    /**
+     * Return a EpetraWrappers::CommunicationPattern and store it for future
      * use.
      */
     EpetraWrappers::CommunicationPattern
index 9a802264816875a5ff8967388abe6c24dc23548f..b4437cfd3fbf9de007806af6aa4048ff3c10a09a 100644 (file)
@@ -34,7 +34,6 @@
 
 #ifdef DEAL_II_WITH_TRILINOS
 #  include <deal.II/lac/trilinos_epetra_communication_pattern.h>
-#  include <deal.II/lac/trilinos_epetra_vector.h>
 #  include <deal.II/lac/trilinos_vector.h>
 
 #  include <Epetra_Import.h>
@@ -507,6 +506,93 @@ namespace LinearAlgebra
 
 
 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::import(
+    const Tpetra::Vector<> &vector,
+    const IndexSet &        source_elements,
+    VectorOperation::values operation,
+    const MPI_Comm &        mpi_comm,
+    const std::shared_ptr<const CommunicationPatternBase>
+      &communication_pattern)
+  {
+    std::shared_ptr<const TpetraWrappers::CommunicationPattern>
+      tpetra_comm_pattern;
+
+    // If no communication pattern is given, create one. Otherwise, use the one
+    // given.
+    if (communication_pattern == nullptr)
+      {
+        // 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.size() == source_stored_elements.size()) &&
+            (source_elements == source_stored_elements))
+          {
+            tpetra_comm_pattern = std::dynamic_pointer_cast<
+              const TpetraWrappers::CommunicationPattern>(comm_pattern);
+            if (tpetra_comm_pattern == nullptr)
+              tpetra_comm_pattern =
+                std::make_shared<const TpetraWrappers::CommunicationPattern>(
+                  create_tpetra_comm_pattern(source_elements, mpi_comm));
+          }
+        else
+          tpetra_comm_pattern =
+            std::make_shared<const TpetraWrappers::CommunicationPattern>(
+              create_tpetra_comm_pattern(source_elements, mpi_comm));
+      }
+    else
+      {
+        tpetra_comm_pattern =
+          std::dynamic_pointer_cast<const TpetraWrappers::CommunicationPattern>(
+            communication_pattern);
+        AssertThrow(tpetra_comm_pattern != nullptr,
+                    ExcMessage(
+                      std::string("The communication pattern is not of type ") +
+                      "LinearAlgebra::TpetraWrappers::CommunicationPattern."));
+      }
+
+    Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export());
+
+    Tpetra::Vector<> target_vector(tpetra_export.getSourceMap());
+    target_vector.doImport(vector, tpetra_export, Tpetra::REPLACE);
+
+    const auto *new_values = target_vector.getData().get();
+    const auto  size       = target_vector.getLocalLength();
+
+    using size_type = std::decay<decltype(size)>::type;
+
+    Assert(size == 0 || values != nullptr, ExcInternalError("Export failed."));
+    AssertDimension(size, stored_elements.n_elements());
+
+    if (operation == VectorOperation::insert)
+      {
+        for (size_type i = 0; i < size; ++i)
+          values[i] = new_values[i];
+      }
+    else if (operation == VectorOperation::add)
+      {
+        for (size_type i = 0; i < size; ++i)
+          values[i] += new_values[i];
+      }
+    else if (operation == VectorOperation::min)
+      {
+        for (size_type i = 0; i < size; ++i)
+          if (std::real(new_values[i]) - std::real(values[i]) < 0.0)
+            values[i] = new_values[i];
+      }
+    else if (operation == VectorOperation::max)
+      {
+        for (size_type i = 0; i < size; ++i)
+          if (std::real(new_values[i]) - std::real(values[i]) > 0.0)
+            values[i] = new_values[i];
+      }
+    else
+      AssertThrow(false, ExcNotImplemented());
+  }
+
+
+
   template <typename Number>
   void
   ReadWriteVector<Number>::import(
@@ -652,6 +738,23 @@ namespace LinearAlgebra
 
 
 
+  template <typename Number>
+  void
+  ReadWriteVector<Number>::import(
+    const LinearAlgebra::TpetraWrappers::Vector &trilinos_vec,
+    VectorOperation::values                      operation,
+    const std::shared_ptr<const CommunicationPatternBase>
+      &communication_pattern)
+  {
+    import(trilinos_vec.trilinos_vector(),
+           trilinos_vec.locally_owned_elements(),
+           operation,
+           trilinos_vec.get_mpi_communicator(),
+           communication_pattern);
+  }
+
+
+
   template <typename Number>
   void
   ReadWriteVector<Number>::import(
@@ -793,6 +896,22 @@ namespace LinearAlgebra
 
 
 #if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
+  template <typename Number>
+  TpetraWrappers::CommunicationPattern
+  ReadWriteVector<Number>::create_tpetra_comm_pattern(
+    const IndexSet &source_index_set,
+    const MPI_Comm &mpi_comm)
+  {
+    source_stored_elements = source_index_set;
+    TpetraWrappers::CommunicationPattern epetra_comm_pattern(
+      source_stored_elements, stored_elements, mpi_comm);
+    comm_pattern = std::make_shared<TpetraWrappers::CommunicationPattern>(
+      source_stored_elements, stored_elements, mpi_comm);
+
+    return epetra_comm_pattern;
+  }
+
+
   template <typename Number>
   EpetraWrappers::CommunicationPattern
   ReadWriteVector<Number>::create_epetra_comm_pattern(
index 25478951f673c370aa0073afa751dc75b9664bcd..006dc31712195ad6f1326f79cd90ef373d39bd8c 100644 (file)
@@ -27,6 +27,7 @@
 #    include <deal.II/lac/exceptions.h>
 #    include <deal.II/lac/full_matrix.h>
 #    include <deal.II/lac/trilinos_epetra_vector.h>
+#    include <deal.II/lac/trilinos_tpetra_vector.h>
 #    include <deal.II/lac/trilinos_vector.h>
 #    include <deal.II/lac/vector_memory.h>
 #    include <deal.II/lac/vector_operation.h>
diff --git a/include/deal.II/lac/trilinos_tpetra_communication_pattern.h b/include/deal.II/lac/trilinos_tpetra_communication_pattern.h
new file mode 100644 (file)
index 0000000..8dc7cf7
--- /dev/null
@@ -0,0 +1,106 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_trilinos_tpetra_communication_pattern_h
+#define dealii_trilinos_tpetra_communication_pattern_h
+
+
+#include <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#  ifdef DEAL_II_WITH_MPI
+
+#    include <deal.II/lac/communication_pattern_base.h>
+
+#    include <Tpetra_Import.hpp>
+
+#    include <memory>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    /**
+     * This class implements a wrapper to Tpetra::Import and Tpetra::Export.
+     */
+    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.
+       */
+      virtual 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.
+       */
+      virtual const MPI_Comm &
+      get_mpi_communicator() const override;
+
+      /**
+       * Return the underlying Tpetra::Import object.
+       */
+      const Tpetra::Import<> &
+      get_tpetra_import() const;
+
+      /**
+       * Return the underlying Tpetra::Export object.
+       */
+      const Tpetra::Export<> &
+      get_tpetra_export() const;
+
+    private:
+      /**
+       * Shared pointer to the MPI communicator used.
+       */
+      std::shared_ptr<const MPI_Comm> comm;
+
+      /**
+       * Shared pointer to the Tpetra::Import object used.
+       */
+      std::unique_ptr<Tpetra::Import<>> tpetra_import;
+
+      /**
+       * Shared pointer to the Tpetra::Export object used.
+       */
+      std::unique_ptr<Tpetra::Export<>> tpetra_export;
+    };
+  } // end of namespace TpetraWrappers
+} // end of namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif
+
+#endif
+
+#endif
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h
new file mode 100644 (file)
index 0000000..df0b0c9
--- /dev/null
@@ -0,0 +1,406 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_trilinos_tpetra_vector_h
+#define dealii_trilinos_tpetra_vector_h
+
+
+#include <deal.II/base/config.h>
+
+#if defined(DEAL_II_WITH_TRILINOS) && defined(DEAL_II_WITH_MPI)
+
+#  include <deal.II/base/index_set.h>
+#  include <deal.II/base/subscriptor.h>
+
+#  include <deal.II/lac/trilinos_tpetra_communication_pattern.h>
+#  include <deal.II/lac/vector_operation.h>
+#  include <deal.II/lac/vector_space_vector.h>
+#  include <deal.II/lac/vector_type_traits.h>
+
+#  include <Teuchos_Comm.hpp>
+#  include <Teuchos_OrdinalTraits.hpp>
+#  include <Tpetra_Core.hpp>
+#  include <Tpetra_Vector.hpp>
+#  include <Tpetra_Version.hpp>
+#  include <mpi.h>
+
+#  include <memory>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  // Forward declaration
+  template <typename Number>
+  class ReadWriteVector;
+
+  namespace TpetraWrappers
+  {
+    /**
+     * This class implements a wrapper to the Trilinos distributed vector
+     * class Tpetra::Vector. This class is derived from the
+     * LinearAlgebra::VectorSpaceVector class. Note however that Tpetra only
+     * works with Number = double.
+     *
+     * @ingroup TrilinosWrappers
+     * @ingroup Vectors
+     * @author Daniel Arndt, 2018
+     */
+    class Vector : public VectorSpaceVector<double>, public Subscriptor
+    {
+    public:
+      /**
+       * Constructor. Create a vector of dimension zero.
+       */
+      Vector();
+
+      /**
+       * Copy constructor. Sets the dimension and the partitioning to that of
+       * the given vector and copies all elements.
+       */
+      Vector(const Vector &V);
+
+      /**
+       * This constructor takes an IndexSet that defines how to distribute the
+       * individual components among the MPI processors. Since it also
+       * includes information about the size of the vector, this is all we
+       * need to generate a %parallel vector.
+       */
+      explicit Vector(const IndexSet &parallel_partitioner,
+                      const MPI_Comm &communicator);
+
+      /**
+       * Reinit functionality. This function destroys the old vector content
+       * and generates a new one based on the input partitioning. The flag
+       * <tt>omit_zeroing_entries</tt> determines whether the vector should be
+       * filled with zero (false) or left untouched (true).
+       */
+      void
+      reinit(const IndexSet &parallel_partitioner,
+             const MPI_Comm &communicator,
+             const bool      omit_zeroing_entries = false);
+
+      /**
+       * Change the dimension to that of the vector V. The elements of V are not
+       * copied.
+       */
+      virtual void
+      reinit(const VectorSpaceVector<double> &V,
+             const bool omit_zeroing_entries = false) override;
+
+      /**
+       * Copy function. This function takes a Vector and copies all the
+       * elements. The Vector will have the same parallel distribution as @p
+       * V.
+       */
+      Vector &
+      operator=(const Vector &V);
+
+      /**
+       * Sets all elements of the vector to the scalar @p s. This operation is
+       * only allowed if @p s is equal to zero.
+       */
+      virtual Vector &
+      operator=(const double s) override;
+
+      /**
+       * 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<double> &                 V,
+        VectorOperation::values                         operation,
+        std::shared_ptr<const CommunicationPatternBase> communication_pattern =
+          std::shared_ptr<const CommunicationPatternBase>()) override;
+
+      /**
+       * Multiply the entire vector by a fixed factor.
+       */
+      virtual Vector &
+      operator*=(const double factor) override;
+
+      /**
+       * Divide the entire vector by a fixed factor.
+       */
+      virtual Vector &
+      operator/=(const double factor) override;
+
+      /**
+       * Add the vector @p V to the present one.
+       */
+      virtual Vector &
+      operator+=(const VectorSpaceVector<double> &V) override;
+
+      /**
+       * Subtract the vector @p V from the present one.
+       */
+      virtual Vector &
+      operator-=(const VectorSpaceVector<double> &V) override;
+
+      /**
+       * Return the scalar product of two vectors. The vectors need to have the
+       * same layout.
+       */
+      virtual double
+      operator*(const VectorSpaceVector<double> &V) const override;
+
+      /**
+       * Add @p a to all components. Note that @p is a scalar not a vector.
+       */
+      virtual void
+      add(const double a) override;
+
+      /**
+       * Simple addition of a multiple of a vector, i.e. <tt>*this +=
+       * a*V</tt>. The vectors need to have the same layout.
+       */
+      virtual void
+      add(const double a, const VectorSpaceVector<double> &V) override;
+
+      /**
+       * Multiple addition of multiple of a vector, i.e. <tt>*this> +=
+       * a*V+b*W</tt>. The vectors need to have the same layout.
+       */
+      virtual void
+      add(const double                     a,
+          const VectorSpaceVector<double> &V,
+          const double                     b,
+          const VectorSpaceVector<double> &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 double                     s,
+           const double                     a,
+           const VectorSpaceVector<double> &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-assignment) by a diagonal scaling matrix. The
+       * vectors need to have the same layout.
+       */
+      virtual void
+      scale(const VectorSpaceVector<double> &scaling_factors) override;
+
+      /**
+       * Assignment <tt>*this = a*V</tt>.
+       */
+      virtual void
+      equ(const double a, const VectorSpaceVector<double> &V) override;
+
+      /**
+       * Return whether the vector contains only elements with value zero.
+       */
+      virtual bool
+      all_zero() const override;
+
+      /**
+       * Return the mean value of the element of this vector.
+       */
+      virtual double
+      mean_value() const override;
+
+      /**
+       * 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 double
+      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 double
+      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 double
+      linfty_norm() const override;
+
+      /**
+       * Performs a combined operation of a vector addition and a subsequent
+       * inner product, returning the value of the inner product. In other
+       * words, the result of this function is the same as if the user called
+       * @code
+       * this->add(a, V);
+       * return_value = *this * W;
+       * @endcode
+       *
+       * The reason this function exists is that this operation involves less
+       * memory transfer than calling the two functions separately. This
+       * method only needs to load three vectors, @p this, @p V, @p W, whereas
+       * calling separate methods means to load the calling vector @p this
+       * twice. Since most vector operations are memory transfer limited, this
+       * reduces the time by 25\% (or 50\% if @p W equals @p this).
+       *
+       * The vectors need to have the same layout.
+       *
+       * For complex-valued vectors, the scalar product in the second step is
+       * implemented as
+       * $\left<v,w\right>=\sum_i v_i \bar{w_i}$.
+       */
+      virtual double
+      add_and_dot(const double                     a,
+                  const VectorSpaceVector<double> &V,
+                  const VectorSpaceVector<double> &W) override;
+      /**
+       * This function always returns false and is present only for backward
+       * compatibility.
+       */
+      bool
+      has_ghost_elements() const;
+
+      /**
+       * Return the global size of the vector, equal to the sum of the number of
+       * locally owned indices among all processors.
+       */
+      virtual size_type
+      size() const override;
+
+      /**
+       * Return the MPI communicator object in use with this object.
+       */
+      MPI_Comm
+      get_mpi_communicator() const;
+
+      /**
+       * Return an index set that describes which elements of this vector are
+       * owned by the current processor. As a consequence, the index sets
+       * returned on different processors if this is a distributed vector will
+       * form disjoint sets that add up to the complete index set. Obviously, if
+       * a vector is created on only one processor, then the result would
+       * satisfy
+       * @code
+       *  vec.locally_owned_elements() == complete_index_set(vec.size())
+       * @endcode
+       */
+      virtual ::dealii::IndexSet
+      locally_owned_elements() const override;
+
+      /**
+       * Return a const reference to the underlying Trilinos
+       * Tpetra::Vector class.
+       */
+      const Tpetra::Vector<> &
+      trilinos_vector() const;
+
+      /**
+       * Return a (modifyable) reference to the underlying Trilinos
+       * Tpetra::Vector class.
+       */
+      Tpetra::Vector<> &
+      trilinos_vector();
+
+      /**
+       * Prints the vector to the output stream @p out.
+       */
+      virtual void
+      print(std::ostream &     out,
+            const unsigned int precision  = 3,
+            const bool         scientific = true,
+            const bool         across     = true) const override;
+
+      /**
+       * Return the memory consumption of this class in bytes.
+       */
+      virtual std::size_t
+      memory_consumption() const override;
+
+      /**
+       * The vectors have different partitioning, i.e. they have use different
+       * IndexSet.
+       */
+      DeclException0(ExcDifferentParallelPartitioning);
+
+      /**
+       * Attempt to perform an operation between two incompatible vector types.
+       *
+       * @ingroup Exceptions
+       */
+      DeclException0(ExcVectorTypeNotCompatible);
+
+      /**
+       * Exception thrown by an error in Trilinos.
+       *
+       * @ingroup Exceptions
+       */
+      DeclException1(ExcTrilinosError,
+                     int,
+                     << "An error with error number " << arg1
+                     << " occurred while calling a Trilinos function");
+
+    private:
+      /**
+       * Create the CommunicationPattern for the communication between the
+       * IndexSet @p source_index_set and the current vector based
+       * on the communicator @p mpi_comm.
+       */
+      void
+      create_tpetra_comm_pattern(const IndexSet &source_index_set,
+                                 const MPI_Comm &mpi_comm);
+
+      /**
+       * Pointer to the actual Tpetra vector object.
+       */
+      std::unique_ptr<Tpetra::Vector<>> 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::shared_ptr<const TpetraWrappers::CommunicationPattern>
+        tpetra_comm_pattern;
+    };
+
+
+    inline bool
+    Vector::has_ghost_elements() const
+    {
+      return false;
+    }
+  } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+
+
+/**
+ * Declare dealii::LinearAlgebra::TpetraWrappers::Vector as distributed vector.
+ */
+template <>
+struct is_serial_vector<LinearAlgebra::TpetraWrappers::Vector> : std::false_type
+{};
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif
+
+#endif
index 791f6894c375cb21b6d591f981d9e4a016a9fa1b..b168260588851e9152285588e41de38e0d70c416 100644 (file)
@@ -18,6 +18,7 @@
 
 
 #include <deal.II/lac/trilinos_epetra_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -108,7 +109,6 @@ namespace internal
     vector[0][trilinos_i] = value;
   }
 
-
   template <>
   inline double
   ElementAccess<LinearAlgebra::EpetraWrappers::Vector>::get(
@@ -122,6 +122,74 @@ namespace internal
 
     return vector[0][trilinos_i];
   }
+
+
+
+  template <>
+  inline void
+  ElementAccess<LinearAlgebra::TpetraWrappers::Vector>::add(
+    const double                           value,
+    const types::global_dof_index          i,
+    LinearAlgebra::TpetraWrappers::Vector &V)
+  {
+    // Extract local indices in the vector.
+    Tpetra::Vector<>                  vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.getMap()->getLocalElement(
+        static_cast<TrilinosWrappers::types::int_type>(i));
+
+    vector.sync<Kokkos::HostSpace>();
+    auto vector_2d = vector.getLocalView<Kokkos::HostSpace>();
+    auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+    // We're going to modify the data on host.
+    vector.modify<Kokkos::HostSpace>();
+    vector_1d(trilinos_i) += value;
+    vector.sync<Tpetra::Vector<>::device_type::memory_space>();
+  }
+
+
+
+  template <>
+  inline void
+  ElementAccess<LinearAlgebra::TpetraWrappers::Vector>::set(
+    const double                           value,
+    const types::global_dof_index          i,
+    LinearAlgebra::TpetraWrappers::Vector &V)
+  {
+    // Extract local indices in the vector.
+    Tpetra::Vector<>                  vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.getMap()->getLocalElement(
+        static_cast<TrilinosWrappers::types::int_type>(i));
+
+    vector.sync<Kokkos::HostSpace>();
+    auto vector_2d = vector.getLocalView<Kokkos::HostSpace>();
+    auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+    // We're going to modify the data on host.
+    vector.modify<Kokkos::HostSpace>();
+    vector_1d(trilinos_i) = value;
+    vector.sync<Tpetra::Vector<>::device_type::memory_space>();
+  }
+
+
+  template <>
+  inline double
+  ElementAccess<LinearAlgebra::TpetraWrappers::Vector>::get(
+    const LinearAlgebra::TpetraWrappers::Vector &V,
+    const types::global_dof_index                i)
+  {
+    // Extract local indices in the vector.
+    Tpetra::Vector<>                  vector = V.trilinos_vector();
+    TrilinosWrappers::types::int_type trilinos_i =
+      vector.getMap()->getLocalElement(
+        static_cast<TrilinosWrappers::types::int_type>(i));
+
+    vector.sync<Kokkos::HostSpace>();
+    auto vector_2d = vector.getLocalView<Kokkos::HostSpace>();
+    auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+    // We're going to modify the data on host.
+    return vector_1d(trilinos_i);
+  }
 #endif
 } // namespace internal
 
index a88bf028846b9a7f9d3abcfd047a3dbf9e945a59..e4b4f3db2ac911ac493d75fca5a7bf760750e45b 100644 (file)
@@ -132,6 +132,38 @@ namespace internal
   };
 
 #  ifdef DEAL_II_WITH_MPI
+  template <>
+  struct MatrixSelector<dealii::LinearAlgebra::TpetraWrappers::Vector>
+  {
+    using Sparsity = ::dealii::TrilinosWrappers::SparsityPattern;
+    using Matrix   = ::dealii::TrilinosWrappers::SparseMatrix;
+
+    static const bool requires_distributed_sparsity_pattern = false;
+
+    template <typename SparsityPatternType, typename DoFHandlerType>
+    static void
+    reinit(Matrix &matrix,
+           Sparsity &,
+           int                        level,
+           const SparsityPatternType &sp,
+           DoFHandlerType &           dh)
+    {
+      const parallel::Triangulation<DoFHandlerType::dimension,
+                                    DoFHandlerType::space_dimension>
+        *dist_tria = dynamic_cast<
+          const parallel::Triangulation<DoFHandlerType::dimension,
+                                        DoFHandlerType::space_dimension> *>(
+          &(dh.get_triangulation()));
+      MPI_Comm communicator =
+        dist_tria != nullptr ? dist_tria->get_communicator() : MPI_COMM_SELF;
+      matrix.reinit(dh.locally_owned_mg_dofs(level + 1),
+                    dh.locally_owned_mg_dofs(level),
+                    sp,
+                    communicator,
+                    true);
+    }
+  };
+
   template <>
   struct MatrixSelector<dealii::LinearAlgebra::EpetraWrappers::Vector>
   {
index 30e5fa5f1fe5e531e2d97c05a731c8a973ac0d1b..aebdafb33768869b55e601a2aae267386c35f874 100644 (file)
@@ -70,6 +70,7 @@
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_memory.h>
index 3e7f0914cf35542914fea0eb0c7c8642a1f732a7..6fe2920ffe390fc3e098d513b517e229d0ad62cd 100644 (file)
@@ -520,6 +520,72 @@ IndexSet::fill_index_vector(std::vector<size_type> &indices) const
 
 #ifdef DEAL_II_WITH_TRILINOS
 
+Tpetra::Map<>
+IndexSet::make_tpetra_map(const MPI_Comm &communicator,
+                          const bool      overlapping) const
+{
+  compress();
+  (void)communicator;
+
+#  ifdef DEBUG
+  if (!overlapping)
+    {
+      const size_type n_global_elements =
+        Utilities::MPI::sum(n_elements(), communicator);
+      Assert(n_global_elements == size(),
+             ExcMessage("You are trying to create an Tpetra::Map object "
+                        "that partitions elements of an index set "
+                        "between processors. However, the union of the "
+                        "index sets on different processors does not "
+                        "contain all indices exactly once: the sum of "
+                        "the number of entries the various processors "
+                        "want to store locally is " +
+                        Utilities::to_string(n_global_elements) +
+                        " whereas the total size of the object to be "
+                        "allocated is " +
+                        Utilities::to_string(size()) +
+                        ". In other words, there are "
+                        "either indices that are not spoken for "
+                        "by any processor, or there are indices that are "
+                        "claimed by multiple processors."));
+    }
+#  endif
+
+  // Find out if the IndexSet is ascending and 1:1. This corresponds to a
+  // linear Tpetra::Map. Overlapping IndexSets are never 1:1.
+  const bool linear =
+    overlapping ? false : is_ascending_and_one_to_one(communicator);
+  if (linear)
+    return Tpetra::Map<>(size(),
+                         n_elements(),
+                         0,
+#  ifdef DEAL_II_WITH_MPI
+                         Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
+#  else
+                         Teuchos::rcp(new Teuchos::Comm<int>())
+#  endif
+    );
+  else
+    {
+      std::vector<size_type> indices;
+      fill_index_vector(indices);
+      std::vector<int> int_indices(indices.size());
+      std::copy(indices.begin(), indices.end(), int_indices.begin());
+      return Tpetra::Map<>(size(),
+                           (n_elements() > 0 ? int_indices.data() : nullptr),
+                           n_elements(),
+                           0,
+#  ifdef DEAL_II_WITH_MPI
+                           Teuchos::rcp(new Teuchos::MpiComm<int>(communicator))
+#  else
+                           Teuchos::rcp(new Teuchos::Comm<int>())
+#  endif
+      );
+    }
+}
+
+
+
 Epetra_Map
 IndexSet::make_trilinos_map(const MPI_Comm &communicator,
                             const bool      overlapping) const
index 9bb02320eb6317ae2f46fe6352fc14b355e29542..3e4d6bed678bf3a73af83ab1b8e88bcec0df0361 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
index 27b798cc124b618695f124f4d6fef235a01276bd..ff4d0d9573224027df4a67dc5e32826d9846fb79 100644 (file)
@@ -979,6 +979,20 @@ namespace Utilities
 
 
 
+    const Teuchos::RCP<const Teuchos::Comm<int>> &
+    tpetra_comm_self()
+    {
+#  ifdef DEAL_II_WITH_MPI
+      static auto communicator = Teuchos::RCP<const Teuchos::Comm<int>>(
+        new Teuchos::MpiComm<int>(MPI_COMM_SELF));
+#  else
+      static auto communicator =
+        Teuchos::RCP<const Teuchos::Comm<int>>(new Teuchos::Comm<int>());
+#  endif
+
+      return communicator;
+    }
+
     const Epetra_Comm &
     comm_self()
     {
index 4f1bdc2b03162341b6bd3a6e8be0fb687457a585..53b47d2238adbf5fb58cc05caf2a4f5105e5569a 100644 (file)
@@ -32,6 +32,7 @@
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
index 556dcd48c8f756408cccad8f49c1d2ae26aa8117..a116ab0858e9cc5eee34fe283177f49cdee7e828 100644 (file)
@@ -32,6 +32,7 @@
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/sparse_matrix.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
index 920f57a61ae041aedc61ea79f15e7f6b6c4c3bee..b2c44908a6b89774e5f5e1c8d38d608d8315c3d5 100644 (file)
@@ -39,6 +39,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_element_access.h>
index 5289c9d1e7ce1bc0a17758f10a57fb3bfea5be35..1659ccfb3a657acbce7f7cd29f6f077989386f2a 100644 (file)
@@ -44,6 +44,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
index a6ce60360fe800f43e20058fbea5b33e3a567aad..973837219e1ba091e958f74259a7d48b9b7c9e88 100644 (file)
@@ -132,6 +132,8 @@ IF(DEAL_II_WITH_TRILINOS)
     trilinos_solver.cc
     trilinos_sparse_matrix.cc
     trilinos_sparsity_pattern.cc
+    trilinos_tpetra_communication_pattern.cc
+    trilinos_tpetra_vector.cc
     trilinos_vector.cc
   )
 
index 6d432ae184b8be2086068ff55c3f3d3502f35201..7bf9e30c3e869f4602321b817e245d041d43e19a 100644 (file)
@@ -22,6 +22,7 @@
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/solver.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_memory.h>
index 291929c5c359c52d35ef95bf140999a4fe6eb75c..350006da3ccabd5811807a7a3b7ebff91ac59961 100644 (file)
@@ -103,6 +103,36 @@ namespace TrilinosWrappers
     {
       return V.trilinos_vector()[0] + V.trilinos_vector().MyLength();
     }
+
+    template <>
+    double *
+    begin(LinearAlgebra::TpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector().getDataNonConst().get();
+    }
+
+    template <>
+    const double *
+    begin(const LinearAlgebra::TpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector().getData().get();
+    }
+
+    template <>
+    double *
+    end(LinearAlgebra::TpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector().getDataNonConst().get() +
+             V.trilinos_vector().getLocalLength();
+    }
+
+    template <>
+    const double *
+    end(const LinearAlgebra::TpetraWrappers::Vector &V)
+    {
+      return V.trilinos_vector().getData().get() +
+             V.trilinos_vector().getLocalLength();
+    }
 #  endif
   } // namespace internal
 
@@ -3454,6 +3484,11 @@ namespace TrilinosWrappers
     dealii::LinearAlgebra::distributed::Vector<double> &,
     const dealii::LinearAlgebra::distributed::Vector<double> &) const;
 #  ifdef DEAL_II_WITH_MPI
+  template void
+  SparseMatrix::vmult(
+    dealii::LinearAlgebra::TpetraWrappers::Vector &,
+    const dealii::LinearAlgebra::TpetraWrappers::Vector &) const;
+
   template void
   SparseMatrix::vmult(
     dealii::LinearAlgebra::EpetraWrappers::Vector &,
@@ -3469,6 +3504,11 @@ namespace TrilinosWrappers
     dealii::LinearAlgebra::distributed::Vector<double> &,
     const dealii::LinearAlgebra::distributed::Vector<double> &) const;
 #  ifdef DEAL_II_WITH_MPI
+  template void
+  SparseMatrix::Tvmult(
+    dealii::LinearAlgebra::TpetraWrappers::Vector &,
+    const dealii::LinearAlgebra::TpetraWrappers::Vector &) const;
+
   template void
   SparseMatrix::Tvmult(
     dealii::LinearAlgebra::EpetraWrappers::Vector &,
@@ -3484,6 +3524,11 @@ namespace TrilinosWrappers
     dealii::LinearAlgebra::distributed::Vector<double> &,
     const dealii::LinearAlgebra::distributed::Vector<double> &) const;
 #  ifdef DEAL_II_WITH_MPI
+  template void
+  SparseMatrix::vmult_add(
+    dealii::LinearAlgebra::TpetraWrappers::Vector &,
+    const dealii::LinearAlgebra::TpetraWrappers::Vector &) const;
+
   template void
   SparseMatrix::vmult_add(
     dealii::LinearAlgebra::EpetraWrappers::Vector &,
@@ -3499,6 +3544,11 @@ namespace TrilinosWrappers
     dealii::LinearAlgebra::distributed::Vector<double> &,
     const dealii::LinearAlgebra::distributed::Vector<double> &) const;
 #  ifdef DEAL_II_WITH_MPI
+  template void
+  SparseMatrix::Tvmult_add(
+    dealii::LinearAlgebra::TpetraWrappers::Vector &,
+    const dealii::LinearAlgebra::TpetraWrappers::Vector &) const;
+
   template void
   SparseMatrix::Tvmult_add(
     dealii::LinearAlgebra::EpetraWrappers::Vector &,
diff --git a/source/lac/trilinos_tpetra_communication_pattern.cc b/source/lac/trilinos_tpetra_communication_pattern.cc
new file mode 100644 (file)
index 0000000..e78dfda
--- /dev/null
@@ -0,0 +1,104 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/std_cxx14/memory.h>
+
+#include <deal.II/lac/trilinos_tpetra_communication_pattern.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#  ifdef DEAL_II_WITH_MPI
+
+#    include <deal.II/base/index_set.h>
+
+#    include <Tpetra_Map.hpp>
+
+#    include <memory>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    CommunicationPattern::CommunicationPattern(
+      const IndexSet &vector_space_vector_index_set,
+      const IndexSet &read_write_vector_index_set,
+      const MPI_Comm &communicator)
+    {
+      // virtual functions called in constructors and destructors never use the
+      // override in a derived class
+      // for clarity be explicit on which function is called
+      CommunicationPattern::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);
+
+      auto vector_space_vector_map = Teuchos::rcp(new Tpetra::Map<>(
+        vector_space_vector_index_set.make_tpetra_map(*comm, false)));
+      auto read_write_vector_map   = Teuchos::rcp(new Tpetra::Map<>(
+        read_write_vector_index_set.make_tpetra_map(*comm, true)));
+
+      // Target map is read_write_vector_map
+      // Source map is vector_space_vector_map. This map must have uniquely
+      // owned GID.
+      tpetra_import =
+        std_cxx14::make_unique<Tpetra::Import<>>(read_write_vector_map,
+                                                 vector_space_vector_map);
+      tpetra_export =
+        std_cxx14::make_unique<Tpetra::Export<>>(read_write_vector_map,
+                                                 vector_space_vector_map);
+    }
+
+
+
+    const MPI_Comm &
+    CommunicationPattern::get_mpi_communicator() const
+    {
+      return *comm;
+    }
+
+
+
+    const Tpetra::Import<> &
+    CommunicationPattern::get_tpetra_import() const
+    {
+      return *tpetra_import;
+    }
+
+
+
+    const Tpetra::Export<> &
+    CommunicationPattern::get_tpetra_export() const
+    {
+      return *tpetra_export;
+    }
+  } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif
+
+#endif
diff --git a/source/lac/trilinos_tpetra_vector.cc b/source/lac/trilinos_tpetra_vector.cc
new file mode 100644 (file)
index 0000000..0a52aa1
--- /dev/null
@@ -0,0 +1,632 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/base/std_cxx14/memory.h>
+
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+#  ifdef DEAL_II_WITH_MPI
+
+#    include <deal.II/base/index_set.h>
+
+#    include <deal.II/lac/read_write_vector.h>
+
+#    include <boost/io/ios_state.hpp>
+
+#    include <Teuchos_DefaultMpiComm.hpp>
+#    include <Tpetra_Import_decl.hpp>
+#    include <Tpetra_Map_decl.hpp>
+
+#    include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    Vector::Vector()
+      : vector(new Tpetra::Vector<>(Teuchos::RCP<Tpetra::Map<>>(
+          new Tpetra::Map<>(0, 0, Utilities::Trilinos::tpetra_comm_self()))))
+    {}
+
+
+
+    Vector::Vector(const Vector &V)
+      : Subscriptor()
+      , vector(new Tpetra::Vector<>(V.trilinos_vector(), Teuchos::Copy))
+    {}
+
+
+
+    Vector::Vector(const IndexSet &parallel_partitioner,
+                   const MPI_Comm &communicator)
+      : vector(new Tpetra::Vector<>(Teuchos::rcp(new Tpetra::Map<>(
+          parallel_partitioner.make_tpetra_map(communicator, false)))))
+    {}
+
+
+
+    void
+    Vector::reinit(const IndexSet &parallel_partitioner,
+                   const MPI_Comm &communicator,
+                   const bool      omit_zeroing_entries)
+    {
+      Tpetra::Map<> input_map =
+        parallel_partitioner.make_tpetra_map(communicator, false);
+      if (vector->getMap()->isSameAs(input_map) == false)
+        vector = std_cxx14::make_unique<Tpetra::Vector<>>(
+          Teuchos::rcp(new Tpetra::Map<>(input_map)));
+      else if (omit_zeroing_entries == false)
+        {
+          vector->putScalar(0.);
+        }
+    }
+
+
+
+    void
+    Vector::reinit(const VectorSpaceVector<double> &V,
+                   const bool                       omit_zeroing_entries)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+
+      reinit(down_V.locally_owned_elements(),
+             down_V.get_mpi_communicator(),
+             omit_zeroing_entries);
+    }
+
+
+
+    Vector &
+    Vector::operator=(const Vector &V)
+    {
+      // Distinguish three cases:
+      //  - First case: both vectors have the same layout.
+      //  - Second case: both vectors have the same size but different layout.
+      //  - Third case: the vectors have different size.
+      if (vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())))
+        *vector = V.trilinos_vector();
+      else
+        {
+          if (size() == V.size())
+            {
+              Tpetra::Import<> data_exchange(vector->getMap(),
+                                             V.trilinos_vector().getMap());
+
+              vector->doImport(V.trilinos_vector(),
+                               data_exchange,
+                               Tpetra::REPLACE);
+            }
+          else
+            vector =
+              std_cxx14::make_unique<Tpetra::Vector<>>(V.trilinos_vector());
+        }
+
+      return *this;
+    }
+
+
+
+    Vector &
+    Vector::operator=(const double s)
+    {
+      Assert(s == 0., ExcMessage("Only 0 can be assigned to a vector."));
+
+      vector->putScalar(s);
+
+      return *this;
+    }
+
+
+
+    void
+    Vector::import(
+      const ReadWriteVector<double> &                 V,
+      VectorOperation::values                         operation,
+      std::shared_ptr<const CommunicationPatternBase> communication_pattern)
+    {
+      // If no communication pattern is given, create one. Otherwsie, use the
+      // one given.
+      if (communication_pattern == nullptr)
+        {
+          // 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.size() !=
+               V.get_stored_elements().size()) ||
+              (source_stored_elements != V.get_stored_elements()))
+            {
+              const Teuchos::MpiComm<int> *mpi_comm =
+                dynamic_cast<const Teuchos::MpiComm<int> *>(
+                  vector->getMap()->getComm().get());
+              Assert(mpi_comm != nullptr, ExcInternalError());
+              create_tpetra_comm_pattern(V.get_stored_elements(),
+                                         *(mpi_comm->getRawMpiComm())());
+            }
+        }
+      else
+        {
+          tpetra_comm_pattern = std::dynamic_pointer_cast<
+            const TpetraWrappers::CommunicationPattern>(communication_pattern);
+          AssertThrow(
+            tpetra_comm_pattern != nullptr,
+            ExcMessage(
+              std::string("The communication pattern is not of type ") +
+              "LinearAlgebra::TpetraWrappers::CommunicationPattern."));
+        }
+
+      Tpetra::Export<> tpetra_export(tpetra_comm_pattern->get_tpetra_export());
+      Tpetra::Vector<> source_vector(tpetra_export.getSourceMap());
+
+      source_vector.sync<Kokkos::HostSpace>();
+      auto x_2d = source_vector.getLocalView<Kokkos::HostSpace>();
+      auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0);
+      source_vector.modify<Kokkos::HostSpace>();
+      const size_t localLength = source_vector.getLocalLength();
+      auto         values_it   = V.begin();
+      for (size_t k = 0; k < localLength; ++k)
+        x_1d(k) = *values_it++;
+      source_vector.sync<Tpetra::Vector<double>::device_type::memory_space>();
+      if (operation == VectorOperation::insert)
+        vector->doExport(source_vector, tpetra_export, Tpetra::REPLACE);
+      else if (operation == VectorOperation::add)
+        vector->doExport(source_vector, tpetra_export, Tpetra::ADD);
+      else
+        AssertThrow(false, ExcNotImplemented());
+    }
+
+
+
+    Vector &
+    Vector::operator*=(const double factor)
+    {
+      AssertIsFinite(factor);
+      vector->scale(factor);
+
+      return *this;
+    }
+
+
+
+    Vector &
+    Vector::operator/=(const double factor)
+    {
+      AssertIsFinite(factor);
+      Assert(factor != 0., ExcZero());
+      *this *= 1. / factor;
+
+      return *this;
+    }
+
+
+
+    Vector &
+    Vector::operator+=(const VectorSpaceVector<double> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      // If the maps are the same we can Update right away.
+      if (vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())))
+        {
+          vector->update(1., down_V.trilinos_vector(), 1.);
+        }
+      else
+        {
+          Assert(this->size() == down_V.size(),
+                 ExcDimensionMismatch(this->size(), down_V.size()));
+
+          // TODO: The code doesn't work as expected so we use a workaround.
+          /*Tpetra::Export<> data_exchange(vector->getMap(),
+                                         down_V.trilinos_vector().getMap());
+          vector->doExport(down_V.trilinos_vector(),
+                           data_exchange,
+                           Tpetra::ADD);*/
+
+          Tpetra::Vector<> dummy(vector->getMap(), false);
+          Tpetra::Import<> data_exchange(dummy.getMap(),
+                                         down_V.trilinos_vector().getMap());
+
+          dummy.doExport(down_V.trilinos_vector(),
+                         data_exchange,
+                         Tpetra::REPLACE);
+
+          vector->update(1.0, dummy, 1.0);
+        }
+
+      return *this;
+    }
+
+
+
+    Vector &
+    Vector::operator-=(const VectorSpaceVector<double> &V)
+    {
+      this->add(-1., V);
+
+      return *this;
+    }
+
+
+
+    double Vector::operator*(const VectorSpaceVector<double> &V) const
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      Assert(this->size() == down_V.size(),
+             ExcDimensionMismatch(this->size(), down_V.size()));
+      Assert(vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()),
+             ExcDifferentParallelPartitioning());
+
+      return vector->dot(down_V.trilinos_vector());
+    }
+
+
+
+    void
+    Vector::add(const double a)
+    {
+      AssertIsFinite(a);
+
+      vector->sync<Kokkos::HostSpace>();
+      auto vector_2d = vector->getLocalView<Kokkos::HostSpace>();
+      auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+      vector->modify<Kokkos::HostSpace>();
+      const size_t localLength = vector->getLocalLength();
+      for (size_t k = 0; k < localLength; ++k)
+        {
+          vector_1d(k) += a;
+        }
+      vector->sync<Tpetra::Vector<double>::device_type::memory_space>();
+    }
+
+
+
+    void
+    Vector::add(const double a, const VectorSpaceVector<double> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      AssertIsFinite(a);
+      Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+
+      vector->update(a, down_V.trilinos_vector(), 1.);
+    }
+
+
+
+    void
+    Vector::add(const double                     a,
+                const VectorSpaceVector<double> &V,
+                const double                     b,
+                const VectorSpaceVector<double> &W)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&W) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      // Downcast W. If fails, throws an exception.
+      const Vector &down_W = dynamic_cast<const Vector &>(W);
+      Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+      Assert(vector->getMap()->isSameAs(*(down_W.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+      AssertIsFinite(a);
+      AssertIsFinite(b);
+
+      vector->update(
+        a, down_V.trilinos_vector(), b, down_W.trilinos_vector(), 1.);
+    }
+
+
+
+    void
+    Vector::sadd(const double                     s,
+                 const double                     a,
+                 const VectorSpaceVector<double> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      *this *= s;
+      // Downcast V. It fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      Vector        tmp(down_V);
+      tmp *= a;
+      *this += tmp;
+    }
+
+
+
+    void
+    Vector::scale(const VectorSpaceVector<double> &scaling_factors)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&scaling_factors) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast scaling_factors. If fails, throws an exception.
+      const Vector &down_scaling_factors =
+        dynamic_cast<const Vector &>(scaling_factors);
+      Assert(vector->getMap()->isSameAs(
+               *(down_scaling_factors.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+
+      vector->elementWiseMultiply(1.,
+                                  *down_scaling_factors.vector,
+                                  *vector,
+                                  0.);
+    }
+
+
+
+    void
+    Vector::equ(const double a, const VectorSpaceVector<double> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector &down_V = dynamic_cast<const Vector &>(V);
+      // If we don't have the same map, copy.
+      if (vector->getMap()->isSameAs(*down_V.trilinos_vector().getMap()) ==
+          false)
+        this->sadd(0., a, V);
+      else
+        {
+          // Otherwise, just update
+          vector->update(a, down_V.trilinos_vector(), 0.);
+        }
+    }
+
+
+
+    bool
+    Vector::all_zero() const
+    {
+      // get a representation of the vector and
+      // loop over all the elements
+      double *      start_ptr = vector->getDataNonConst().get();
+      const double *ptr       = start_ptr,
+                   *eptr      = start_ptr + vector->getLocalLength();
+      unsigned int flag       = 0;
+      while (ptr != eptr)
+        {
+          if (*ptr != 0)
+            {
+              flag = 1;
+              break;
+            }
+          ++ptr;
+        }
+
+      // Check that the vector is zero on _all_ processors.
+      const Teuchos::MpiComm<int> *mpi_comm =
+        dynamic_cast<const Teuchos::MpiComm<int> *>(
+          vector->getMap()->getComm().get());
+      Assert(mpi_comm != nullptr, ExcInternalError());
+      unsigned int num_nonzero =
+        Utilities::MPI::sum(flag, *(mpi_comm->getRawMpiComm())());
+
+      return num_nonzero == 0;
+    }
+
+
+
+    double
+    Vector::mean_value() const
+    {
+      return vector->meanValue();
+    }
+
+
+
+    double
+    Vector::l1_norm() const
+    {
+      return vector->norm1();
+    }
+
+
+
+    double
+    Vector::l2_norm() const
+    {
+      return vector->norm2();
+    }
+
+
+
+    double
+    Vector::linfty_norm() const
+    {
+      return vector->normInf();
+    }
+
+
+
+    double
+    Vector::add_and_dot(const double                     a,
+                        const VectorSpaceVector<double> &V,
+                        const VectorSpaceVector<double> &W)
+    {
+      this->add(a, V);
+
+      return *this * W;
+    }
+
+
+
+    Vector::size_type
+    Vector::size() const
+    {
+      return vector->getGlobalLength();
+    }
+
+
+
+    MPI_Comm
+    Vector::get_mpi_communicator() const
+    {
+      const auto tpetra_comm = dynamic_cast<const Teuchos::MpiComm<int> *>(
+        vector->getMap()->getComm().get());
+      Assert(tpetra_comm != nullptr, ExcInternalError());
+      return *(tpetra_comm->getRawMpiComm())();
+    }
+
+
+
+    ::dealii::IndexSet
+    Vector::locally_owned_elements() const
+    {
+      IndexSet is(size());
+
+      // easy case: local range is contiguous
+      if (vector->getMap()->isContiguous())
+        {
+#    ifndef DEAL_II_WITH_64BIT_INDICES
+          is.add_range(vector->getMap()->getMinGlobalIndex(),
+                       vector->getMap()->getMaxGlobalIndex() + 1);
+#    else
+          is.add_range(vector->getMap()->getMinGlobalIndex(),
+                       vector->getMap()->getMaxGlobalIndex() + 1);
+#    endif
+        }
+      else if (vector->getLocalLength() > 0)
+        {
+          const size_type n_indices = vector->getLocalLength();
+          auto vector_indices       = vector->getMap()->getMyGlobalIndices();
+          is.add_indices((unsigned int *)&vector_indices[0],
+                         (unsigned int *)&vector_indices[0] + n_indices);
+        }
+      is.compress();
+
+      return is;
+    }
+
+
+
+    const Tpetra::Vector<> &
+    Vector::trilinos_vector() const
+    {
+      return *vector;
+    }
+
+
+
+    Tpetra::Vector<> &
+    Vector::trilinos_vector()
+    {
+      return *vector;
+    }
+
+
+
+    void
+    Vector::print(std::ostream &     out,
+                  const unsigned int precision,
+                  const bool         scientific,
+                  const bool         across) const
+    {
+      AssertThrow(out, ExcIO());
+      boost::io::ios_flags_saver restore_flags(out);
+
+      // Get a representation of the vector and loop over all
+      // the elements
+      const auto val = vector->get1dView();
+
+      out.precision(precision);
+      if (scientific)
+        out.setf(std::ios::scientific, std::ios::floatfield);
+      else
+        out.setf(std::ios::fixed, std::ios::floatfield);
+
+      vector->sync<Kokkos::HostSpace>();
+      auto         vector_2d    = vector->getLocalView<Kokkos::HostSpace>();
+      auto         vector_1d    = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+      const size_t local_length = vector->getLocalLength();
+
+      if (across)
+        for (unsigned int i = 0; i < local_length; ++i)
+          out << i << "->" << vector->getMap()->getGlobalElement(i) << ": "
+              << vector_1d(i) << std::endl;
+      //          out << val[i] << ' ';
+      else
+        for (unsigned int i = 0; i < local_length; ++i)
+          out << val[i] << std::endl;
+      out << std::endl;
+
+      // restore the representation
+      // of the vector
+      AssertThrow(out, ExcIO());
+    }
+
+
+
+    std::size_t
+    Vector::memory_consumption() const
+    {
+      return sizeof(*this) +
+             vector->getLocalLength() *
+               (sizeof(double) + sizeof(TrilinosWrappers::types::int_type));
+    }
+
+
+
+    void
+    Vector::create_tpetra_comm_pattern(const IndexSet &source_index_set,
+                                       const MPI_Comm &mpi_comm)
+    {
+      source_stored_elements = source_index_set;
+      tpetra_comm_pattern =
+        std::make_shared<TpetraWrappers::CommunicationPattern>(
+          locally_owned_elements(), source_index_set, mpi_comm);
+    }
+  } // namespace TpetraWrappers
+} // namespace LinearAlgebra
+
+DEAL_II_NAMESPACE_CLOSE
+
+#  endif
+
+#endif
index f15f32800e75473214b10af2fba3a1785b5b9f4e..3fa5affd9bd93a0bcec185fec15d8fe6a94aaf48 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_memory.templates.h>
index 8744c6cf40499804b08e3f6a4663ebf4bb3f947a..a32b68a096889eae076e4ba885de730793e64c85 100644 (file)
@@ -19,6 +19,7 @@
 #include <deal.II/lac/la_vector.h>
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
index e085c3d3825be4fcdebaefef73da53dd485dd4c3..676a33714ecbde11f4f0654556c2282a7db5f128 100644 (file)
@@ -21,6 +21,7 @@
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_epetra_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 #include <deal.II/lac/vector_element_access.h>
index 73f574bdc116ad04df31146b2a199ec0244095a6..d062b4ec2bd9d0dcebdeb2f98b39abb6dcf5b171 100644 (file)
@@ -39,6 +39,7 @@
 #include <deal.II/lac/petsc_block_vector.h>
 #include <deal.II/lac/petsc_vector.h>
 #include <deal.II/lac/trilinos_parallel_block_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
 #include <deal.II/lac/trilinos_vector.h>
 #include <deal.II/lac/vector.h>
 
diff --git a/tests/trilinos/tpetra_vector_01.cc b/tests/trilinos/tpetra_vector_01.cc
new file mode 100644 (file)
index 0000000..c10ac87
--- /dev/null
@@ -0,0 +1,225 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+// Check LinearAlgebra::TpetraWrappers::Vector assignment and import
+
+
+void
+test()
+{
+  IndexSet     parallel_partitioner_1(10);
+  IndexSet     parallel_partitioner_2(10);
+  unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  if (rank == 0)
+    {
+      parallel_partitioner_1.add_range(0, 5);
+      parallel_partitioner_2.add_range(0, 3);
+    }
+  else
+    {
+      parallel_partitioner_1.add_range(5, 10);
+      parallel_partitioner_2.add_range(3, 10);
+    }
+  parallel_partitioner_1.compress();
+  parallel_partitioner_2.compress();
+  LinearAlgebra::TpetraWrappers::Vector a;
+  LinearAlgebra::TpetraWrappers::Vector b(parallel_partitioner_1,
+                                          MPI_COMM_WORLD);
+  LinearAlgebra::TpetraWrappers::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, MPI_COMM_WORLD);
+  AssertThrow(a.size() == 10, ExcMessage("Vector has the wrong size."));
+
+  AssertThrow(parallel_partitioner_1 == b.locally_owned_elements(),
+              ExcMessage("IndexSet has been modified."));
+  AssertThrow(parallel_partitioner_2 == a.locally_owned_elements(),
+              ExcMessage("IndexSet has been modified."));
+
+  IndexSet read_write_index_set(10);
+  if (rank == 0)
+    read_write_index_set.add_range(0, 6);
+  else
+    read_write_index_set.add_range(4, 10);
+  read_write_index_set.compress();
+
+  LinearAlgebra::ReadWriteVector<double> read_write_1(read_write_index_set);
+  LinearAlgebra::ReadWriteVector<double> read_write_2(read_write_index_set);
+  LinearAlgebra::ReadWriteVector<double> read_write_3(read_write_index_set);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        {
+          read_write_1[i] = i;
+          read_write_2[i] = 5. + i;
+        }
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        {
+          read_write_1[i] = i;
+          read_write_2[i] = 5. + i;
+        }
+    }
+
+  a.import(read_write_2, VectorOperation::insert);
+  AssertThrow(a.size() == 10, ExcMessage("Vector has the wrong size."));
+
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        {
+          AssertThrow(read_write_2[i] == read_write_3[i],
+                      ExcMessage("Vector a has been modified."));
+        }
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(read_write_2[i] == read_write_3[i],
+                    ExcMessage("Vector a has been modified."));
+    }
+
+  b.import(read_write_1, VectorOperation::insert);
+  AssertThrow(b.size() == 10, ExcMessage("Vector has the wrong size."));
+
+  read_write_3.import(b, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(read_write_1[i] == read_write_3[i],
+                    ExcMessage("Vector b has been modified."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(read_write_1[i] == read_write_3[i],
+                    ExcMessage("Vector b has been modified."));
+    }
+
+  c.import(read_write_2, VectorOperation::insert);
+  AssertThrow(c.size() == 10, ExcMessage("Vector has the wrong size."));
+
+  read_write_3.import(c, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(read_write_2[i] == read_write_3[i],
+                    ExcMessage("Vector c has been modified."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(read_write_2[i] == read_write_3[i],
+                    ExcMessage("Vector c has been modified."));
+    }
+
+  a *= 2;
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(2. * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in operator *=."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(2. * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in operator *=."));
+    }
+
+  c /= 2.;
+  read_write_3.import(c, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(0.5 * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in operator /=."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(0.5 * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in operator /=."));
+    }
+
+  b += a;
+  read_write_3.import(b, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(2. * read_write_2[i] + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in operator +=."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(2. * read_write_2[i] + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in operator +=."));
+    }
+
+  b -= c;
+  read_write_3.import(b, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 6; ++i)
+        AssertThrow(1.5 * read_write_2[i] + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in operator -=."));
+    }
+  else
+    {
+      for (unsigned int i = 4; i < 10; ++i)
+        AssertThrow(1.5 * read_write_2[i] + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in operator -=."));
+    }
+
+  b.import(read_write_1, VectorOperation::insert);
+  c.import(read_write_1, VectorOperation::insert);
+  const double val = b * c;
+  AssertThrow(val == 285., ExcMessage("Problem in operator *."));
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  test();
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/trilinos/tpetra_vector_01.mpirun=2.output b/tests/trilinos/tpetra_vector_01.mpirun=2.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/tpetra_vector_02.cc b/tests/trilinos/tpetra_vector_02.cc
new file mode 100644 (file)
index 0000000..8e48753
--- /dev/null
@@ -0,0 +1,216 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 - 2017 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+// Check LinearAlgebra::TpetraWrappers::Vector add and sadd.
+
+void
+test()
+{
+  IndexSet     parallel_partitioner_1(10);
+  IndexSet     parallel_partitioner_2(10);
+  unsigned int rank = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  if (rank == 0)
+    {
+      parallel_partitioner_1.add_range(0, 5);
+      parallel_partitioner_2.add_range(0, 3);
+    }
+  else
+    {
+      parallel_partitioner_1.add_range(5, 10);
+      parallel_partitioner_2.add_range(3, 10);
+    }
+  parallel_partitioner_1.compress();
+  parallel_partitioner_2.compress();
+  LinearAlgebra::TpetraWrappers::Vector a(parallel_partitioner_1,
+                                          MPI_COMM_WORLD);
+  LinearAlgebra::TpetraWrappers::Vector b(parallel_partitioner_1,
+                                          MPI_COMM_WORLD);
+  LinearAlgebra::TpetraWrappers::Vector c(parallel_partitioner_2,
+                                          MPI_COMM_WORLD);
+
+  IndexSet read_write_index_set(10);
+  if (rank == 0)
+    read_write_index_set.add_range(0, 5);
+  else
+    read_write_index_set.add_range(5, 10);
+  read_write_index_set.compress();
+
+  LinearAlgebra::ReadWriteVector<double> read_write_1(read_write_index_set);
+  LinearAlgebra::ReadWriteVector<double> read_write_2(read_write_index_set);
+  LinearAlgebra::ReadWriteVector<double> read_write_3(read_write_index_set);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        {
+          read_write_1[i] = i;
+          read_write_2[i] = 5. + i;
+        }
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        {
+          read_write_1[i] = i;
+          read_write_2[i] = 5. + i;
+        }
+    }
+
+  a.import(read_write_1, VectorOperation::insert);
+  b.import(read_write_2, VectorOperation::insert);
+  c.import(read_write_2, VectorOperation::insert);
+
+  a.add(1.);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(1. + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in add(scalar)."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(1. + read_write_1[i] == read_write_3[i],
+                    ExcMessage("Problem in add(scalar)."));
+    }
+
+  a.add(2., b);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(1. + read_write_1[i] + 2. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in add(scalar,Vector)."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(1. + read_write_1[i] + 2. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in add(scalar,Vector)."));
+    }
+
+
+  LinearAlgebra::TpetraWrappers::Vector d(a);
+  a.add(2., b, 3., d);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(4. + 4. * read_write_1[i] + 10. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(4. + 4. * read_write_1[i] + 10. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in add(scalar,Vector,scalar,Vector)."));
+    }
+
+
+  a.import(read_write_1, VectorOperation::insert);
+  a.sadd(3., 2., c);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(3. * read_write_1[i] + 2. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(3. * read_write_1[i] + 2. * read_write_2[i] ==
+                      read_write_3[i],
+                    ExcMessage("Problem in sadd(scalar,scalar,Vector)."));
+    }
+
+
+  a.import(read_write_1, VectorOperation::insert);
+  a.scale(b);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(read_write_1[i] * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in scale."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(read_write_1[i] * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in scale."));
+    }
+
+
+  a.equ(2., c);
+  read_write_3.import(a, VectorOperation::insert);
+  if (rank == 0)
+    {
+      for (unsigned int i = 0; i < 5; ++i)
+        AssertThrow(2. * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in scale."));
+    }
+  else
+    {
+      for (unsigned int i = 5; i < 10; ++i)
+        AssertThrow(2. * read_write_2[i] == read_write_3[i],
+                    ExcMessage("Problem in equ."));
+    }
+
+
+  AssertThrow(b.l1_norm() == 95., ExcMessage("Problem in l1_norm."));
+
+  const double eps = 1e-6;
+  AssertThrow(std::fabs(b.l2_norm() - 31.3847096) < eps,
+              ExcMessage("Problem in l2_norm"));
+
+  AssertThrow(b.linfty_norm() == 14., ExcMessage("Problem in linfty_norm."));
+
+  a.import(read_write_1, VectorOperation::insert);
+  const double val = a.add_and_dot(2., a, b);
+  AssertThrow(val == 1530., ExcMessage("Problem in add_and_dot"));
+}
+
+
+int
+main(int argc, char **argv)
+{
+  initlog();
+  deallog.depth_console(0);
+
+  Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1);
+
+  test();
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/trilinos/tpetra_vector_02.mpirun=2.output b/tests/trilinos/tpetra_vector_02.mpirun=2.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK
diff --git a/tests/trilinos/tpetra_vector_03.cc b/tests/trilinos/tpetra_vector_03.cc
new file mode 100644 (file)
index 0000000..4927aae
--- /dev/null
@@ -0,0 +1,101 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 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.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+// Check LinearAlgebra::TpetraWrappers::Vector::import
+// for VectorOperation::add and check LinearAlgebra::ReadWriteVector
+// for VectorOperation::add/min/max.
+
+#include <deal.II/base/index_set.h>
+#include <deal.II/base/utilities.h>
+
+#include <deal.II/lac/read_write_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#include <iostream>
+#include <vector>
+
+#include "../tests.h"
+
+void
+test()
+{
+  unsigned int my_id   = Utilities::MPI::this_mpi_process(MPI_COMM_WORLD);
+  unsigned int n_procs = Utilities::MPI::n_mpi_processes(MPI_COMM_WORLD);
+
+  IndexSet locally_owned(n_procs * 2);
+  locally_owned.add_range(my_id * 2, my_id * 2 + 2);
+  locally_owned.compress();
+
+  LinearAlgebra::TpetraWrappers::Vector v(locally_owned, MPI_COMM_WORLD);
+
+  IndexSet     workaround_set(n_procs * 2);
+  unsigned int my_first_index  = (my_id * 2 + 2) % (n_procs * 2);
+  unsigned int my_second_index = my_id * 2 + 1;
+  workaround_set.add_index(my_first_index);
+  workaround_set.add_index(my_second_index);
+  workaround_set.add_index(0);
+  workaround_set.compress();
+  LinearAlgebra::ReadWriteVector<double> rw_vector(workaround_set);
+
+  rw_vector(my_first_index)  = my_id + 10;
+  rw_vector(my_second_index) = my_id + 100;
+  rw_vector(0)               = 1.;
+  // rw_vector(2) = 1.;
+  v.import(rw_vector, VectorOperation::add);
+  deallog << "Tpetra first import add:" << std::endl;
+  v.print(deallog.get_file_stream());
+  rw_vector.print(deallog.get_file_stream());
+
+  rw_vector(my_first_index)  = my_id + 20;
+  rw_vector(my_second_index) = my_id + 200;
+  rw_vector(0)               = 2.;
+  // rw_vector(2) = 3.;
+  v.import(rw_vector, VectorOperation::add);
+  deallog << "Tpetra second import add:" << std::endl;
+  v.print(deallog.get_file_stream());
+  rw_vector.print(deallog.get_file_stream());
+
+  rw_vector.import(v, VectorOperation::add);
+  deallog << "ReadWrite import add:" << std::endl;
+  rw_vector.print(deallog.get_file_stream());
+
+  rw_vector(my_first_index)  = my_id + 100;
+  rw_vector(my_second_index) = 1;
+  rw_vector(0)               = 4.;
+  // rw_vector(2) = 3.;
+  rw_vector.import(v, VectorOperation::min);
+  deallog << "ReadWrite import min:" << std::endl;
+  rw_vector.print(deallog.get_file_stream());
+
+  rw_vector(my_first_index)  = my_id + 100;
+  rw_vector(my_second_index) = 1;
+  rw_vector(0)               = 4.;
+  // rw_vector(2) = 3.;
+  rw_vector.import(v, VectorOperation::max);
+  deallog << "ReadWrite import max:" << std::endl;
+  rw_vector.print(deallog.get_file_stream());
+}
+
+
+
+int
+main(int argc, char **argv)
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(
+    argc, argv, testing_max_num_threads());
+
+  MPILogInitAll log;
+  test();
+}
diff --git a/tests/trilinos/tpetra_vector_03.mpirun=2.output b/tests/trilinos/tpetra_vector_03.mpirun=2.output
new file mode 100644 (file)
index 0000000..5d5f9c0
--- /dev/null
@@ -0,0 +1,70 @@
+
+DEAL:0::Tpetra first import add:
+0->0: 2.000e+00
+1->1: 1.000e+02
+
+IndexSet: {[0,2]}
+
+[0]: 1.000e+00
+[1]: 1.000e+02
+[2]: 1.000e+01
+DEAL:0::Tpetra second import add:
+0->0: 4.000e+00
+1->1: 2.000e+02
+
+IndexSet: {[0,2]}
+
+[0]: 2.000e+00
+[1]: 2.000e+02
+[2]: 2.000e+01
+DEAL:0::ReadWrite import add:
+IndexSet: {[0,2]}
+
+[0]: 6.000e+00
+[1]: 4.000e+02
+[2]: 5.000e+01
+DEAL:0::ReadWrite import min:
+IndexSet: {[0,2]}
+
+[0]: 4.000e+00
+[1]: 1.000e+00
+[2]: 3.000e+01
+DEAL:0::ReadWrite import max:
+IndexSet: {[0,2]}
+
+[0]: 4.000e+00
+[1]: 2.000e+02
+[2]: 1.000e+02
+
+DEAL:1::Tpetra first import add:
+0->2: 1.000e+01
+1->3: 1.010e+02
+
+IndexSet: {0, 3}
+
+[0]: 1.000e+00
+[3]: 1.010e+02
+DEAL:1::Tpetra second import add:
+0->2: 3.000e+01
+1->3: 2.010e+02
+
+IndexSet: {0, 3}
+
+[0]: 2.000e+00
+[3]: 2.010e+02
+DEAL:1::ReadWrite import add:
+IndexSet: {0, 3}
+
+[0]: 6.000e+00
+[3]: 4.020e+02
+DEAL:1::ReadWrite import min:
+IndexSet: {0, 3}
+
+[0]: 4.000e+00
+[3]: 1.000e+00
+DEAL:1::ReadWrite import max:
+IndexSet: {0, 3}
+
+[0]: 4.000e+00
+[3]: 2.010e+02
+
diff --git a/tests/trilinos/tpetra_vector_03.mpirun=4.output b/tests/trilinos/tpetra_vector_03.mpirun=4.output
new file mode 100644 (file)
index 0000000..bb8d36d
--- /dev/null
@@ -0,0 +1,146 @@
+
+DEAL:0::Tpetra first import add:
+0->0: 4.000e+00
+1->1: 1.000e+02
+
+IndexSet: {[0,2]}
+
+[0]: 1.000e+00
+[1]: 1.000e+02
+[2]: 1.000e+01
+DEAL:0::Tpetra second import add:
+0->0: 8.000e+00
+1->1: 2.000e+02
+
+IndexSet: {[0,2]}
+
+[0]: 2.000e+00
+[1]: 2.000e+02
+[2]: 2.000e+01
+DEAL:0::ReadWrite import add:
+IndexSet: {[0,2]}
+
+[0]: 1.000e+01
+[1]: 4.000e+02
+[2]: 5.000e+01
+DEAL:0::ReadWrite import min:
+IndexSet: {[0,2]}
+
+[0]: 4.000e+00
+[1]: 1.000e+00
+[2]: 3.000e+01
+DEAL:0::ReadWrite import max:
+IndexSet: {[0,2]}
+
+[0]: 8.000e+00
+[1]: 2.000e+02
+[2]: 1.000e+02
+
+DEAL:1::Tpetra first import add:
+0->2: 1.000e+01
+1->3: 1.010e+02
+
+IndexSet: {0, [3,4]}
+
+[0]: 1.000e+00
+[3]: 1.010e+02
+[4]: 1.100e+01
+DEAL:1::Tpetra second import add:
+0->2: 3.000e+01
+1->3: 2.010e+02
+
+IndexSet: {0, [3,4]}
+
+[0]: 2.000e+00
+[3]: 2.010e+02
+[4]: 2.100e+01
+DEAL:1::ReadWrite import add:
+IndexSet: {0, [3,4]}
+
+[0]: 1.000e+01
+[3]: 4.020e+02
+[4]: 5.300e+01
+DEAL:1::ReadWrite import min:
+IndexSet: {0, [3,4]}
+
+[0]: 4.000e+00
+[3]: 1.000e+00
+[4]: 3.200e+01
+DEAL:1::ReadWrite import max:
+IndexSet: {0, [3,4]}
+
+[0]: 8.000e+00
+[3]: 2.010e+02
+[4]: 1.010e+02
+
+
+DEAL:2::Tpetra first import add:
+0->4: 1.100e+01
+1->5: 1.020e+02
+
+IndexSet: {0, [5,6]}
+
+[0]: 1.000e+00
+[5]: 1.020e+02
+[6]: 1.200e+01
+DEAL:2::Tpetra second import add:
+0->4: 3.200e+01
+1->5: 2.020e+02
+
+IndexSet: {0, [5,6]}
+
+[0]: 2.000e+00
+[5]: 2.020e+02
+[6]: 2.200e+01
+DEAL:2::ReadWrite import add:
+IndexSet: {0, [5,6]}
+
+[0]: 1.000e+01
+[5]: 4.040e+02
+[6]: 5.600e+01
+DEAL:2::ReadWrite import min:
+IndexSet: {0, [5,6]}
+
+[0]: 4.000e+00
+[5]: 1.000e+00
+[6]: 3.400e+01
+DEAL:2::ReadWrite import max:
+IndexSet: {0, [5,6]}
+
+[0]: 8.000e+00
+[5]: 2.020e+02
+[6]: 1.020e+02
+
+
+DEAL:3::Tpetra first import add:
+0->6: 1.200e+01
+1->7: 1.030e+02
+
+IndexSet: {0, 7}
+
+[0]: 1.000e+00
+[7]: 1.030e+02
+DEAL:3::Tpetra second import add:
+0->6: 3.400e+01
+1->7: 2.030e+02
+
+IndexSet: {0, 7}
+
+[0]: 2.000e+00
+[7]: 2.030e+02
+DEAL:3::ReadWrite import add:
+IndexSet: {0, 7}
+
+[0]: 1.000e+01
+[7]: 4.060e+02
+DEAL:3::ReadWrite import min:
+IndexSet: {0, 7}
+
+[0]: 4.000e+00
+[7]: 1.000e+00
+DEAL:3::ReadWrite import max:
+IndexSet: {0, 7}
+
+[0]: 8.000e+00
+[7]: 2.030e+02
+

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.