]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Move implementation to templates file
authorDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 22 Feb 2019 14:04:18 +0000 (15:04 +0100)
committerDaniel Arndt <daniel.arndt@iwr.uni-heidelberg.de>
Fri, 22 Feb 2019 14:04:18 +0000 (15:04 +0100)
include/deal.II/lac/trilinos_tpetra_vector.templates.h [new file with mode: 0644]
source/lac/trilinos_tpetra_vector.cc

diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
new file mode 100644 (file)
index 0000000..b102775
--- /dev/null
@@ -0,0 +1,681 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2018 - 2019 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_templates_h
+#define dealii__trilinos_tpetra_vector_templates_h
+
+#include <deal.II/base/std_cxx14/memory.h>
+
+#include <deal.II/lac/trilinos_tpetra_vector.h>
+
+#ifdef DEAL_II_TRILINOS_WITH_TPETRA
+
+#  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_def.hpp>
+#    include <Tpetra_Map_def.hpp>
+
+#    include <memory>
+
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace LinearAlgebra
+{
+  namespace TpetraWrappers
+  {
+    template <typename Number>
+    Vector<Number>::Vector()
+      : Subscriptor()
+      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
+          Teuchos::RCP<Tpetra::Map<int, types::global_dof_index>>(
+            new Tpetra::Map<int, types::global_dof_index>(
+              0,
+              0,
+              Utilities::Trilinos::tpetra_comm_self()))))
+    {}
+
+
+
+    template <typename Number>
+    Vector<Number>::Vector(const Vector<Number> &V)
+      : Subscriptor()
+      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
+          V.trilinos_vector(),
+          Teuchos::Copy))
+    {}
+
+
+
+    template <typename Number>
+    Vector<Number>::Vector(const IndexSet &parallel_partitioner,
+                           const MPI_Comm &communicator)
+      : Subscriptor()
+      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
+          Teuchos::rcp(new Tpetra::Map<int, types::global_dof_index>(
+            parallel_partitioner.make_tpetra_map(communicator, false)))))
+    {}
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::reinit(const IndexSet &parallel_partitioner,
+                           const MPI_Comm &communicator,
+                           const bool      omit_zeroing_entries)
+    {
+      Tpetra::Map<int, types::global_dof_index> input_map =
+        parallel_partitioner.make_tpetra_map(communicator, false);
+      if (vector->getMap()->isSameAs(input_map) == false)
+        vector = std_cxx14::make_unique<
+          Tpetra::Vector<Number, int, types::global_dof_index>>(Teuchos::rcp(
+          new Tpetra::Map<int, types::global_dof_index>(input_map)));
+      else if (omit_zeroing_entries == false)
+        {
+          vector->putScalar(0.);
+        }
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::reinit(const VectorSpaceVector<Number> &V,
+                           const bool omit_zeroing_entries)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
+
+      reinit(down_V.locally_owned_elements(),
+             down_V.get_mpi_communicator(),
+             omit_zeroing_entries);
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator=(const Vector<Number> &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<int, types::global_dof_index> 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<Number, int, types::global_dof_index>>(
+              V.trilinos_vector());
+        }
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator=(const Number s)
+    {
+      Assert(s == Number(0.),
+             ExcMessage("Only 0 can be assigned to a vector."));
+
+      vector->putScalar(s);
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::import(
+      const ReadWriteVector<Number> &                 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<int, types::global_dof_index> tpetra_export(
+        tpetra_comm_pattern->get_tpetra_export());
+      Tpetra::Vector<Number, int, types::global_dof_index> source_vector(
+        tpetra_export.getSourceMap());
+
+      source_vector.template sync<Kokkos::HostSpace>();
+      auto x_2d = source_vector.template getLocalView<Kokkos::HostSpace>();
+      auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0);
+      source_vector.template 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.template sync<
+        typename Tpetra::Vector<Number, int, types::global_dof_index>::
+          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());
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator*=(const Number factor)
+    {
+      AssertIsFinite(factor);
+      vector->scale(factor);
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator/=(const Number factor)
+    {
+      AssertIsFinite(factor);
+      Assert(factor != Number(0.), ExcZero());
+      *this *= Number(1.) / factor;
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator+=(const VectorSpaceVector<Number> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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: Tpetra doesn't have a combine mode that also updates local
+          // elements, maybe there is a better workaround.
+          Tpetra::Vector<Number, int, types::global_dof_index> dummy(
+            vector->getMap(), false);
+          Tpetra::Import<int, types::global_dof_index> data_exchange(
+            down_V.trilinos_vector().getMap(), dummy.getMap());
+
+          dummy.doImport(down_V.trilinos_vector(),
+                         data_exchange,
+                         Tpetra::INSERT);
+
+          vector->update(1.0, dummy, 1.0);
+        }
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    Vector<Number> &
+    Vector<Number>::operator-=(const VectorSpaceVector<Number> &V)
+    {
+      this->add(-1., V);
+
+      return *this;
+    }
+
+
+
+    template <typename Number>
+    Number Vector<Number>::operator*(const VectorSpaceVector<Number> &V) const
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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());
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::add(const Number a)
+    {
+      AssertIsFinite(a);
+
+      vector->template sync<Kokkos::HostSpace>();
+      auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>();
+      auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
+      vector->template modify<Kokkos::HostSpace>();
+      const size_t localLength = vector->getLocalLength();
+      for (size_t k = 0; k < localLength; ++k)
+        {
+          vector_1d(k) += a;
+        }
+      vector->template sync<
+        typename Tpetra::Vector<Number, int, types::global_dof_index>::
+          device_type::memory_space>();
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
+      AssertIsFinite(a);
+      Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+
+      vector->update(a, down_V.trilinos_vector(), 1.);
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::add(const Number                     a,
+                        const VectorSpaceVector<Number> &V,
+                        const Number                     b,
+                        const VectorSpaceVector<Number> &W)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
+      // Downcast W. If fails, throws an exception.
+      const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(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.);
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::sadd(const Number                     s,
+                         const Number                     a,
+                         const VectorSpaceVector<Number> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      *this *= s;
+      // Downcast V. It fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
+      Vector<Number>        tmp(down_V);
+      tmp *= a;
+      *this += tmp;
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&scaling_factors) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast scaling_factors. If fails, throws an exception.
+      const Vector<Number> &down_scaling_factors =
+        dynamic_cast<const Vector<Number> &>(scaling_factors);
+      Assert(vector->getMap()->isSameAs(
+               *(down_scaling_factors.trilinos_vector().getMap())),
+             ExcDifferentParallelPartitioning());
+
+      vector->elementWiseMultiply(1.,
+                                  *down_scaling_factors.vector,
+                                  *vector,
+                                  0.);
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
+    {
+      // Check that casting will work.
+      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
+             ExcVectorTypeNotCompatible());
+
+      // Downcast V. If fails, throws an exception.
+      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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.);
+        }
+    }
+
+
+
+    template <typename Number>
+    bool
+    Vector<Number>::all_zero() const
+    {
+      // get a representation of the vector and
+      // loop over all the elements
+      Number *      start_ptr = vector->getDataNonConst().get();
+      const Number *ptr       = start_ptr,
+                   *eptr      = start_ptr + vector->getLocalLength();
+      unsigned int flag       = 0;
+      while (ptr != eptr)
+        {
+          if (*ptr != Number(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;
+    }
+
+
+
+    template <typename Number>
+    Number
+    Vector<Number>::mean_value() const
+    {
+      return vector->meanValue();
+    }
+
+
+
+    template <typename Number>
+    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
+    Vector<Number>::l1_norm() const
+    {
+      return vector->norm1();
+    }
+
+
+
+    template <typename Number>
+    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
+    Vector<Number>::l2_norm() const
+    {
+      return vector->norm2();
+    }
+
+
+
+    template <typename Number>
+    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
+    Vector<Number>::linfty_norm() const
+    {
+      return vector->normInf();
+    }
+
+
+
+    template <typename Number>
+    Number
+    Vector<Number>::add_and_dot(const Number                     a,
+                                const VectorSpaceVector<Number> &V,
+                                const VectorSpaceVector<Number> &W)
+    {
+      this->add(a, V);
+
+      return *this * W;
+    }
+
+
+
+    template <typename Number>
+    typename Vector<Number>::size_type
+    Vector<Number>::size() const
+    {
+      return vector->getGlobalLength();
+    }
+
+
+
+    template <typename Number>
+    MPI_Comm
+    Vector<Number>::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())();
+    }
+
+
+
+    template <typename Number>
+    ::dealii::IndexSet
+    Vector<Number>::locally_owned_elements() const
+    {
+      IndexSet is(size());
+
+      // easy case: local range is contiguous
+      if (vector->getMap()->isContiguous())
+        {
+          is.add_range(vector->getMap()->getMinGlobalIndex(),
+                       vector->getMap()->getMaxGlobalIndex() + 1);
+        }
+      else if (vector->getLocalLength() > 0)
+        {
+          const size_type n_indices = vector->getLocalLength();
+          std::vector<types::global_dof_index> vector_indices;
+          vector_indices.reserve(n_indices);
+          for (unsigned int i = 0; i < n_indices; ++i)
+            vector_indices.push_back(vector->getMap()->getGlobalElement(i));
+
+          is.add_indices(vector_indices.data(),
+                         vector_indices.data() + n_indices);
+        }
+      is.compress();
+
+      return is;
+    }
+
+
+
+    template <typename Number>
+    const Tpetra::Vector<Number, int, types::global_dof_index> &
+    Vector<Number>::trilinos_vector() const
+    {
+      return *vector;
+    }
+
+
+
+    template <typename Number>
+    Tpetra::Vector<Number, int, types::global_dof_index> &
+    Vector<Number>::trilinos_vector()
+    {
+      return *vector;
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::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->template sync<Kokkos::HostSpace>();
+      auto vector_2d = vector->template 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 << vector_1d(i) << ' ';
+      else
+        for (unsigned int i = 0; i < local_length; ++i)
+          out << vector_1d(i) << std::endl;
+      out << std::endl;
+
+      // restore the representation
+      // of the vector
+      AssertThrow(out, ExcIO());
+    }
+
+
+
+    template <typename Number>
+    std::size_t
+    Vector<Number>::memory_consumption() const
+    {
+      return sizeof(*this) +
+             vector->getLocalLength() *
+               (sizeof(Number) + sizeof(TrilinosWrappers::types::int_type));
+    }
+
+
+
+    template <typename Number>
+    void
+    Vector<Number>::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
+
+#endif
index ea218f8309f876be3df65c8f0f40cb394d9c4c77..b269bdabcf88a30e8361c0a1ff31aad67a6506aa 100644 (file)
 //
 // ---------------------------------------------------------------------
 
-#include <deal.II/base/std_cxx14/memory.h>
-
-#include <deal.II/lac/trilinos_tpetra_vector.h>
+#include <deal.II/lac/trilinos_tpetra_vector.templates.h>
 
 #ifdef DEAL_II_TRILINOS_WITH_TPETRA
-
 #  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_def.hpp>
-#    include <Tpetra_Map_def.hpp>
-
-#    include <memory>
-
-
 DEAL_II_NAMESPACE_OPEN
 
 namespace LinearAlgebra
 {
   namespace TpetraWrappers
   {
-    template <typename Number>
-    Vector<Number>::Vector()
-      : Subscriptor()
-      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
-          Teuchos::RCP<Tpetra::Map<int, types::global_dof_index>>(
-            new Tpetra::Map<int, types::global_dof_index>(
-              0,
-              0,
-              Utilities::Trilinos::tpetra_comm_self()))))
-    {}
-
-
-
-    template <typename Number>
-    Vector<Number>::Vector(const Vector<Number> &V)
-      : Subscriptor()
-      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
-          V.trilinos_vector(),
-          Teuchos::Copy))
-    {}
-
-
-
-    template <typename Number>
-    Vector<Number>::Vector(const IndexSet &parallel_partitioner,
-                           const MPI_Comm &communicator)
-      : Subscriptor()
-      , vector(new Tpetra::Vector<Number, int, types::global_dof_index>(
-          Teuchos::rcp(new Tpetra::Map<int, types::global_dof_index>(
-            parallel_partitioner.make_tpetra_map(communicator, false)))))
-    {}
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::reinit(const IndexSet &parallel_partitioner,
-                           const MPI_Comm &communicator,
-                           const bool      omit_zeroing_entries)
-    {
-      Tpetra::Map<int, types::global_dof_index> input_map =
-        parallel_partitioner.make_tpetra_map(communicator, false);
-      if (vector->getMap()->isSameAs(input_map) == false)
-        vector = std_cxx14::make_unique<
-          Tpetra::Vector<Number, int, types::global_dof_index>>(Teuchos::rcp(
-          new Tpetra::Map<int, types::global_dof_index>(input_map)));
-      else if (omit_zeroing_entries == false)
-        {
-          vector->putScalar(0.);
-        }
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::reinit(const VectorSpaceVector<Number> &V,
-                           const bool omit_zeroing_entries)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
-
-      reinit(down_V.locally_owned_elements(),
-             down_V.get_mpi_communicator(),
-             omit_zeroing_entries);
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator=(const Vector<Number> &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<int, types::global_dof_index> 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<Number, int, types::global_dof_index>>(
-              V.trilinos_vector());
-        }
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator=(const Number s)
-    {
-      Assert(s == Number(0.),
-             ExcMessage("Only 0 can be assigned to a vector."));
-
-      vector->putScalar(s);
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::import(
-      const ReadWriteVector<Number> &                 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<int, types::global_dof_index> tpetra_export(
-        tpetra_comm_pattern->get_tpetra_export());
-      Tpetra::Vector<Number, int, types::global_dof_index> source_vector(
-        tpetra_export.getSourceMap());
-
-      source_vector.template sync<Kokkos::HostSpace>();
-      auto x_2d = source_vector.template getLocalView<Kokkos::HostSpace>();
-      auto x_1d = Kokkos::subview(x_2d, Kokkos::ALL(), 0);
-      source_vector.template 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.template sync<
-        typename Tpetra::Vector<Number, int, types::global_dof_index>::
-          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());
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator*=(const Number factor)
-    {
-      AssertIsFinite(factor);
-      vector->scale(factor);
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator/=(const Number factor)
-    {
-      AssertIsFinite(factor);
-      Assert(factor != Number(0.), ExcZero());
-      *this *= Number(1.) / factor;
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator+=(const VectorSpaceVector<Number> &V)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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: Tpetra doesn't have a combine mode that also updates local
-          // elements, maybe there is a better workaround.
-          Tpetra::Vector<Number, int, types::global_dof_index> dummy(
-            vector->getMap(), false);
-          Tpetra::Import<int, types::global_dof_index> data_exchange(
-            down_V.trilinos_vector().getMap(), dummy.getMap());
-
-          dummy.doImport(down_V.trilinos_vector(),
-                         data_exchange,
-                         Tpetra::INSERT);
-
-          vector->update(1.0, dummy, 1.0);
-        }
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator-=(const VectorSpaceVector<Number> &V)
-    {
-      this->add(-1., V);
-
-      return *this;
-    }
-
-
-
-    template <typename Number>
-    Number Vector<Number>::operator*(const VectorSpaceVector<Number> &V) const
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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());
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::add(const Number a)
-    {
-      AssertIsFinite(a);
-
-      vector->template sync<Kokkos::HostSpace>();
-      auto vector_2d = vector->template getLocalView<Kokkos::HostSpace>();
-      auto vector_1d = Kokkos::subview(vector_2d, Kokkos::ALL(), 0);
-      vector->template modify<Kokkos::HostSpace>();
-      const size_t localLength = vector->getLocalLength();
-      for (size_t k = 0; k < localLength; ++k)
-        {
-          vector_1d(k) += a;
-        }
-      vector->template sync<
-        typename Tpetra::Vector<Number, int, types::global_dof_index>::
-          device_type::memory_space>();
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::add(const Number a, const VectorSpaceVector<Number> &V)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
-      AssertIsFinite(a);
-      Assert(vector->getMap()->isSameAs(*(down_V.trilinos_vector().getMap())),
-             ExcDifferentParallelPartitioning());
-
-      vector->update(a, down_V.trilinos_vector(), 1.);
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::add(const Number                     a,
-                        const VectorSpaceVector<Number> &V,
-                        const Number                     b,
-                        const VectorSpaceVector<Number> &W)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&W) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
-      // Downcast W. If fails, throws an exception.
-      const Vector<Number> &down_W = dynamic_cast<const Vector<Number> &>(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.);
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::sadd(const Number                     s,
-                         const Number                     a,
-                         const VectorSpaceVector<Number> &V)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      *this *= s;
-      // Downcast V. It fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(V);
-      Vector<Number>        tmp(down_V);
-      tmp *= a;
-      *this += tmp;
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::scale(const VectorSpaceVector<Number> &scaling_factors)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&scaling_factors) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast scaling_factors. If fails, throws an exception.
-      const Vector<Number> &down_scaling_factors =
-        dynamic_cast<const Vector<Number> &>(scaling_factors);
-      Assert(vector->getMap()->isSameAs(
-               *(down_scaling_factors.trilinos_vector().getMap())),
-             ExcDifferentParallelPartitioning());
-
-      vector->elementWiseMultiply(1.,
-                                  *down_scaling_factors.vector,
-                                  *vector,
-                                  0.);
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::equ(const Number a, const VectorSpaceVector<Number> &V)
-    {
-      // Check that casting will work.
-      Assert(dynamic_cast<const Vector<Number> *>(&V) != nullptr,
-             ExcVectorTypeNotCompatible());
-
-      // Downcast V. If fails, throws an exception.
-      const Vector<Number> &down_V = dynamic_cast<const Vector<Number> &>(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.);
-        }
-    }
-
-
-
-    template <typename Number>
-    bool
-    Vector<Number>::all_zero() const
-    {
-      // get a representation of the vector and
-      // loop over all the elements
-      Number *      start_ptr = vector->getDataNonConst().get();
-      const Number *ptr       = start_ptr,
-                   *eptr      = start_ptr + vector->getLocalLength();
-      unsigned int flag       = 0;
-      while (ptr != eptr)
-        {
-          if (*ptr != Number(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;
-    }
-
-
-
-    template <typename Number>
-    Number
-    Vector<Number>::mean_value() const
-    {
-      return vector->meanValue();
-    }
-
-
-
-    template <typename Number>
-    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
-    Vector<Number>::l1_norm() const
-    {
-      return vector->norm1();
-    }
-
-
-
-    template <typename Number>
-    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
-    Vector<Number>::l2_norm() const
-    {
-      return vector->norm2();
-    }
-
-
-
-    template <typename Number>
-    typename LinearAlgebra::VectorSpaceVector<Number>::real_type
-    Vector<Number>::linfty_norm() const
-    {
-      return vector->normInf();
-    }
-
-
-
-    template <typename Number>
-    Number
-    Vector<Number>::add_and_dot(const Number                     a,
-                                const VectorSpaceVector<Number> &V,
-                                const VectorSpaceVector<Number> &W)
-    {
-      this->add(a, V);
-
-      return *this * W;
-    }
-
-
-
-    template <typename Number>
-    typename Vector<Number>::size_type
-    Vector<Number>::size() const
-    {
-      return vector->getGlobalLength();
-    }
-
-
-
-    template <typename Number>
-    MPI_Comm
-    Vector<Number>::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())();
-    }
-
-
-
-    template <typename Number>
-    ::dealii::IndexSet
-    Vector<Number>::locally_owned_elements() const
-    {
-      IndexSet is(size());
-
-      // easy case: local range is contiguous
-      if (vector->getMap()->isContiguous())
-        {
-          is.add_range(vector->getMap()->getMinGlobalIndex(),
-                       vector->getMap()->getMaxGlobalIndex() + 1);
-        }
-      else if (vector->getLocalLength() > 0)
-        {
-          const size_type n_indices = vector->getLocalLength();
-          std::vector<types::global_dof_index> vector_indices;
-          vector_indices.reserve(n_indices);
-          for (unsigned int i = 0; i < n_indices; ++i)
-            vector_indices.push_back(vector->getMap()->getGlobalElement(i));
-
-          is.add_indices(vector_indices.data(),
-                         vector_indices.data() + n_indices);
-        }
-      is.compress();
-
-      return is;
-    }
-
-
-
-    template <typename Number>
-    const Tpetra::Vector<Number, int, types::global_dof_index> &
-    Vector<Number>::trilinos_vector() const
-    {
-      return *vector;
-    }
-
-
-
-    template <typename Number>
-    Tpetra::Vector<Number, int, types::global_dof_index> &
-    Vector<Number>::trilinos_vector()
-    {
-      return *vector;
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::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->template sync<Kokkos::HostSpace>();
-      auto vector_2d = vector->template 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 << vector_1d(i) << ' ';
-      else
-        for (unsigned int i = 0; i < local_length; ++i)
-          out << vector_1d(i) << std::endl;
-      out << std::endl;
-
-      // restore the representation
-      // of the vector
-      AssertThrow(out, ExcIO());
-    }
-
-
-
-    template <typename Number>
-    std::size_t
-    Vector<Number>::memory_consumption() const
-    {
-      return sizeof(*this) +
-             vector->getLocalLength() *
-               (sizeof(Number) + sizeof(TrilinosWrappers::types::int_type));
-    }
-
-
-
-    template <typename Number>
-    void
-    Vector<Number>::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);
-    }
-
     template class Vector<float>;
     template class Vector<double>;
 #    ifdef DEAL_II_WITH_COMPLEX_VALUES
@@ -679,5 +36,4 @@ namespace LinearAlgebra
 DEAL_II_NAMESPACE_CLOSE
 
 #  endif
-
 #endif

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.