From: Daniel Arndt <arndtd@ornl.gov>
Date: Tue, 30 Jan 2024 20:18:17 +0000 (-0500)
Subject: Template TrilinosWrappers::Vector on memory space
X-Git-Tag: relicensing~89^2~1
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=18b2d1ff1373aede34518dabc91556fda158bc9d;p=dealii.git

Template TrilinosWrappers::Vector on memory space
---

diff --git a/include/deal.II/base/template_constraints.h b/include/deal.II/base/template_constraints.h
index 8f5374df5b..bbbd0685fa 100644
--- a/include/deal.II/base/template_constraints.h
+++ b/include/deal.II/base/template_constraints.h
@@ -668,7 +668,7 @@ namespace LinearAlgebra
 #  ifdef DEAL_II_TRILINOS_WITH_TPETRA
   namespace TpetraWrappers
   {
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     class Vector;
 
     template <typename Number, typename MemorySpace>
diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
index 858115ca85..f9d314c772 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
@@ -464,7 +464,8 @@ namespace LinearAlgebra
        * of the matrix.
        */
       void
-      vmult(Vector<Number> &dst, const Vector<Number> &src) const;
+      vmult(Vector<Number, MemorySpace>       &dst,
+            const Vector<Number, MemorySpace> &src) const;
 
       /*
        * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
@@ -474,7 +475,8 @@ namespace LinearAlgebra
        * Source and destination must not be the same vector.
        */
       void
-      Tvmult(Vector<Number> &dst, const Vector<Number> &src) const;
+      Tvmult(Vector<Number, MemorySpace>       &dst,
+             const Vector<Number, MemorySpace> &src) const;
 
       /**
        * Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
@@ -483,7 +485,8 @@ namespace LinearAlgebra
        * Source and destination must not be the same vector.
        */
       void
-      vmult_add(Vector<Number> &dst, const Vector<Number> &src) const;
+      vmult_add(Vector<Number, MemorySpace>       &dst,
+                const Vector<Number, MemorySpace> &src) const;
 
 
       /**
@@ -494,7 +497,8 @@ namespace LinearAlgebra
        * Source and destination must not be the same vector.
        */
       void
-      Tvmult_add(Vector<Number> &dst, const Vector<Number> &src) const;
+      Tvmult_add(Vector<Number, MemorySpace>       &dst,
+                 const Vector<Number, MemorySpace> &src) const;
       /** @} */
 
       /**
diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
index d13fdac5ff..2a1d74fcad 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -585,8 +585,9 @@ namespace LinearAlgebra
 
     template <typename Number, typename MemorySpace>
     void
-    SparseMatrix<Number, MemorySpace>::vmult(Vector<Number>       &dst,
-                                             const Vector<Number> &src) const
+    SparseMatrix<Number, MemorySpace>::vmult(
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
     {
       Assert(&src != &dst, ExcSourceEqualsDestination());
       Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
@@ -601,8 +602,9 @@ namespace LinearAlgebra
 
     template <typename Number, typename MemorySpace>
     void
-    SparseMatrix<Number, MemorySpace>::Tvmult(Vector<Number>       &dst,
-                                              const Vector<Number> &src) const
+    SparseMatrix<Number, MemorySpace>::Tvmult(
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
     {
       Assert(&src != &dst, ExcSourceEqualsDestination());
       Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
@@ -618,8 +620,8 @@ namespace LinearAlgebra
     template <typename Number, typename MemorySpace>
     void
     SparseMatrix<Number, MemorySpace>::vmult_add(
-      Vector<Number>       &dst,
-      const Vector<Number> &src) const
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
     {
       Assert(&src != &dst, ExcSourceEqualsDestination());
       Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
@@ -639,8 +641,8 @@ namespace LinearAlgebra
     template <typename Number, typename MemorySpace>
     void
     SparseMatrix<Number, MemorySpace>::Tvmult_add(
-      Vector<Number>       &dst,
-      const Vector<Number> &src) const
+      Vector<Number, MemorySpace>       &dst,
+      const Vector<Number, MemorySpace> &src) const
     {
       Assert(&src != &dst, ExcSourceEqualsDestination());
       Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.h b/include/deal.II/lac/trilinos_tpetra_vector.h
index 1411617a8f..feea186481 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.h
@@ -125,7 +125,7 @@ namespace LinearAlgebra
        *
        * @ingroup TpetraWrappers
        */
-      template <typename Number>
+      template <typename Number, typename MemorySpace>
       class VectorReference
       {
       private:
@@ -133,7 +133,8 @@ namespace LinearAlgebra
          * Constructor. It is made private so as to only allow the actual vector
          * class to create it.
          */
-        VectorReference(Vector<Number> &vector, const size_type index);
+        VectorReference(Vector<Number, MemorySpace> &vector,
+                        const size_type              index);
 
       public:
         /**
@@ -209,7 +210,7 @@ namespace LinearAlgebra
         /**
          * Point to the vector we are referencing.
          */
-        Vector<Number> &vector;
+        Vector<Number, MemorySpace> &vector;
 
         /**
          * Index of the referenced element of the vector.
@@ -218,7 +219,7 @@ namespace LinearAlgebra
 
         // Make the vector class a friend, so that it can create objects of the
         // present type.
-        friend class Vector<Number>;
+        friend class Vector<Number, MemorySpace>;
       }; // class VectorReference
 
     } // namespace internal
@@ -232,38 +233,40 @@ namespace LinearAlgebra
      * class Tpetra::Vector. This class requires Trilinos to be
      * compiled with MPI support.
      *
-     * Tpetra uses Kokkos for thread-parallelism and chooses the execution and
-     * memory space automatically depending on Kokkos configuration. The
-     * priority is ranked from highest to lowest:
-     * - GPU backend
-     * - host parallel backend
-     * - Kokkos::Serial
-     *
-     * In case Kokkos was configured with GPU support, this class performs its
-     * actions on the GPU. In particular, there is no need for manually
-     * synchronizing memory between host and @ref GlossDevice "device".
+     * Moreover, this class takes an optional template argument for
+     * the memory space used. By default, all memory is allocated on the CPU.
      *
      * @ingroup TpetraWrappers
      * @ingroup Vectors
      */
-    template <typename Number>
+    template <typename Number, typename MemorySpace = dealii::MemorySpace::Host>
     class Vector : public ReadVector<Number>, public Subscriptor
     {
     public:
       /**
        * Declare some of the standard types used in all containers.
        */
+#  if DEAL_II_TRILINOS_VERSION_GTE(14, 2, 0)
+      using NodeType = Tpetra::KokkosCompat::KokkosDeviceWrapperNode<
+        typename MemorySpace::kokkos_space::execution_space,
+        typename MemorySpace::kokkos_space>;
+#  else
+      using NodeType = Kokkos::Compat::KokkosDeviceWrapperNode<
+        typename MemorySpace::kokkos_space::execution_space,
+        typename MemorySpace::kokkos_space>;
+#  endif
       using value_type = Number;
       using real_type  = typename numbers::NumberTraits<Number>::real_type;
       using size_type  = types::global_dof_index;
-      using reference  = internal::VectorReference<Number>;
-      using MapType = Tpetra::Map<int, dealii::types::signed_global_dof_index>;
-      using VectorType =
-        Tpetra::Vector<Number, int, dealii::types::signed_global_dof_index>;
+      using reference  = internal::VectorReference<Number, MemorySpace>;
+      using MapType =
+        Tpetra::Map<int, dealii::types::signed_global_dof_index, NodeType>;
+      using VectorType = Tpetra::
+        Vector<Number, int, dealii::types::signed_global_dof_index, NodeType>;
       using ExportType =
-        Tpetra::Export<int, dealii::types::signed_global_dof_index>;
+        Tpetra::Export<int, dealii::types::signed_global_dof_index, NodeType>;
       using ImportType =
-        Tpetra::Import<int, dealii::types::signed_global_dof_index>;
+        Tpetra::Import<int, dealii::types::signed_global_dof_index, NodeType>;
 
       /**
        * @name 1: Basic Object-handling
@@ -352,7 +355,8 @@ namespace LinearAlgebra
        * copied.
        */
       void
-      reinit(const Vector<Number> &V, const bool omit_zeroing_entries = false);
+      reinit(const Vector<Number, MemorySpace> &V,
+             const bool                         omit_zeroing_entries = false);
 
       /**
        * Extract a range of elements all at once.
@@ -494,20 +498,20 @@ namespace LinearAlgebra
        * Add the vector @p V to the present one.
        */
       Vector &
-      operator+=(const Vector<Number> &V);
+      operator+=(const Vector<Number, MemorySpace> &V);
 
       /**
        * Subtract the vector @p V from the present one.
        */
       Vector &
-      operator-=(const Vector<Number> &V);
+      operator-=(const Vector<Number, MemorySpace> &V);
 
       /**
        * Return the scalar product of two vectors. The vectors need to have the
        * same layout.
        */
       Number
-      operator*(const Vector<Number> &V) const;
+      operator*(const Vector<Number, MemorySpace> &V) const;
 
       /**
        * Add @p a to all components. Note that @p is a scalar not a vector.
@@ -520,17 +524,17 @@ namespace LinearAlgebra
        * a*V</tt>. The vectors need to have the same layout.
        */
       void
-      add(const Number a, const Vector<Number> &V);
+      add(const Number a, const Vector<Number, MemorySpace> &V);
 
       /**
        * Multiple addition of multiple of a vector, i.e. <tt>*this> +=
        * a*V+b*W</tt>. The vectors need to have the same layout.
        */
       void
-      add(const Number          a,
-          const Vector<Number> &V,
-          const Number          b,
-          const Vector<Number> &W);
+      add(const Number                       a,
+          const Vector<Number, MemorySpace> &V,
+          const Number                       b,
+          const Vector<Number, MemorySpace> &W);
 
       /**
        * A collective add operation: This function adds a whole set of values
@@ -556,7 +560,9 @@ namespace LinearAlgebra
        * = s*(*this)+a*V</tt>.
        */
       void
-      sadd(const Number s, const Number a, const Vector<Number> &V);
+      sadd(const Number                       s,
+           const Number                       a,
+           const Vector<Number, MemorySpace> &V);
 
       /**
        * A collective set operation: instead of setting individual elements of a
@@ -576,13 +582,13 @@ namespace LinearAlgebra
        * vectors need to have the same layout.
        */
       void
-      scale(const Vector<Number> &scaling_factors);
+      scale(const Vector<Number, MemorySpace> &scaling_factors);
 
       /**
        * Assignment <tt>*this = a*V</tt>.
        */
       void
-      equ(const Number a, const Vector<Number> &V);
+      equ(const Number a, const Vector<Number, MemorySpace> &V);
 
       /**
        * Return whether the vector contains only elements with value zero.
@@ -648,9 +654,9 @@ namespace LinearAlgebra
        * $\left<v,w\right>=\sum_i v_i \bar{w_i}$.
        */
       Number
-      add_and_dot(const Number          a,
-                  const Vector<Number> &V,
-                  const Vector<Number> &W);
+      add_and_dot(const Number                       a,
+                  const Vector<Number, MemorySpace> &V,
+                  const Vector<Number, MemorySpace> &W);
 
       /** @} */
 
@@ -903,34 +909,34 @@ namespace LinearAlgebra
         tpetra_comm_pattern;
 
       // Make the reference class a friend.
-      friend class internal::VectorReference<Number>;
+      friend class internal::VectorReference<Number, MemorySpace>;
     };
 
 
     /* ------------------------- Inline functions ---------------------- */
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline bool
-    Vector<Number>::has_ghost_elements() const
+    Vector<Number, MemorySpace>::has_ghost_elements() const
     {
       return has_ghost;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline bool
-    Vector<Number>::is_compressed() const
+    Vector<Number, MemorySpace>::is_compressed() const
     {
       return compressed;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline void
-    Vector<Number>::add(const std::vector<size_type>      &indices,
-                        const std::vector<TrilinosScalar> &values)
+    Vector<Number, MemorySpace>::add(const std::vector<size_type>      &indices,
+                                     const std::vector<TrilinosScalar> &values)
     {
       // if we have ghost values, do not allow
       // writing to this vector at all.
@@ -941,11 +947,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline void
-    Vector<Number>::add(const size_type  n_elements,
-                        const size_type *indices,
-                        const Number    *values)
+    Vector<Number, MemorySpace>::add(const size_type  n_elements,
+                                     const size_type *indices,
+                                     const Number    *values)
     {
 #  if DEAL_II_TRILINOS_VERSION_GTE(13, 2, 0)
       auto vector_2d_local = vector->template getLocalView<Kokkos::HostSpace>(
@@ -1021,11 +1027,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline void
-    Vector<Number>::set(const size_type  n_elements,
-                        const size_type *indices,
-                        const Number    *values)
+    Vector<Number, MemorySpace>::set(const size_type  n_elements,
+                                     const size_type *indices,
+                                     const Number    *values)
     {
       // if we have ghost values, do not allow
       // writing to this vector at all.
@@ -1078,23 +1084,23 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    inline internal::VectorReference<Number>
-    Vector<Number>::operator()(const size_type index)
+    template <typename Number, typename MemorySpace>
+    inline internal::VectorReference<Number, MemorySpace>
+    Vector<Number, MemorySpace>::operator()(const size_type index)
     {
       return internal::VectorReference(*this, index);
     }
 
-    template <typename Number>
-    inline internal::VectorReference<Number>
-    Vector<Number>::operator[](const size_type index)
+    template <typename Number, typename MemorySpace>
+    inline internal::VectorReference<Number, MemorySpace>
+    Vector<Number, MemorySpace>::operator[](const size_type index)
     {
       return operator()(index);
     }
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     inline Number
-    Vector<Number>::operator[](const size_type index) const
+    Vector<Number, MemorySpace>::operator[](const size_type index) const
     {
       return operator()(index);
     }
@@ -1104,18 +1110,20 @@ namespace LinearAlgebra
     // VectorReference
     namespace internal
     {
-      template <typename Number>
-      inline VectorReference<Number>::VectorReference(Vector<Number> &vector,
-                                                      const size_type index)
+      template <typename Number, typename MemorySpace>
+      inline VectorReference<Number, MemorySpace>::VectorReference(
+        Vector<Number, MemorySpace> &vector,
+        const size_type              index)
         : vector(vector)
         , index(index)
       {}
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator=(const VectorReference<Number> &r) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator=(
+        const VectorReference<Number, MemorySpace> &r) const
       {
         // as explained in the class
         // documentation, this is not the copy
@@ -1128,9 +1136,10 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline VectorReference<Number> &
-      VectorReference<Number>::operator=(const VectorReference<Number> &r)
+      template <typename Number, typename MemorySpace>
+      inline VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator=(
+        const VectorReference<Number, MemorySpace> &r)
       {
         // as above
         *this = static_cast<Number>(r);
@@ -1140,9 +1149,9 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator=(const Number &value) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator=(const Number &value) const
       {
         vector.set(1, &index, &value);
         return *this;
@@ -1150,9 +1159,10 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator+=(const Number &value) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator+=(
+        const Number &value) const
       {
         vector.add(1, &index, &value);
         return *this;
@@ -1160,9 +1170,10 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator-=(const Number &value) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator-=(
+        const Number &value) const
       {
         Number new_value = -value;
         vector.add(1, &index, &new_value);
@@ -1171,9 +1182,10 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator*=(const Number &value) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator*=(
+        const Number &value) const
       {
         Number new_value = static_cast<Number>(*this) * value;
         vector.set(1, &index, &new_value);
@@ -1182,9 +1194,10 @@ namespace LinearAlgebra
 
 
 
-      template <typename Number>
-      inline const VectorReference<Number> &
-      VectorReference<Number>::operator/=(const Number &value) const
+      template <typename Number, typename MemorySpace>
+      inline const VectorReference<Number, MemorySpace> &
+      VectorReference<Number, MemorySpace>::operator/=(
+        const Number &value) const
       {
         Number new_value = static_cast<Number>(*this) / value;
         vector.set(1, &index, &new_value);
@@ -1203,9 +1216,9 @@ namespace LinearAlgebra
 /**
  * Declare dealii::LinearAlgebra::TpetraWrappers::Vector as distributed vector.
  */
-template <typename Number>
-struct is_serial_vector<LinearAlgebra::TpetraWrappers::Vector<Number>>
-  : std::false_type
+template <typename Number, typename MemorySpace>
+struct is_serial_vector<
+  LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>> : std::false_type
 {};
 
 DEAL_II_NAMESPACE_CLOSE
diff --git a/include/deal.II/lac/trilinos_tpetra_vector.templates.h b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
index ba226f7fbe..8a61db0c0f 100644
--- a/include/deal.II/lac/trilinos_tpetra_vector.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_vector.templates.h
@@ -44,8 +44,8 @@ namespace LinearAlgebra
 {
   namespace TpetraWrappers
   {
-    template <typename Number>
-    Vector<Number>::Vector()
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace>::Vector()
       : Subscriptor()
       , compressed(true)
       , has_ghost(false)
@@ -58,8 +58,8 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number>::Vector(const Vector<Number> &V)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace>::Vector(const Vector<Number, MemorySpace> &V)
       : Subscriptor()
       , compressed(V.compressed)
       , has_ghost(V.has_ghost)
@@ -69,8 +69,8 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number>::Vector(const Teuchos::RCP<VectorType> V)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace>::Vector(const Teuchos::RCP<VectorType> V)
       : Subscriptor()
       , compressed(true)
       , has_ghost(V->getMap()->isOneToOne() == false)
@@ -79,9 +79,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number>::Vector(const IndexSet &parallel_partitioner,
-                           const MPI_Comm  communicator)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace>::Vector(const IndexSet &parallel_partitioner,
+                                        const MPI_Comm  communicator)
       : Subscriptor()
       , compressed(true)
       , has_ghost(false)
@@ -91,11 +91,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number>::Vector(const IndexSet &locally_owned_entries,
-                           const IndexSet &ghost_entries,
-                           const MPI_Comm  communicator,
-                           const bool      vector_writable)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace>::Vector(const IndexSet &locally_owned_entries,
+                                        const IndexSet &ghost_entries,
+                                        const MPI_Comm  communicator,
+                                        const bool      vector_writable)
       : Subscriptor()
     {
       if (!vector_writable)
@@ -135,11 +135,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::reinit(const IndexSet &parallel_partitioner,
-                           const MPI_Comm  communicator,
-                           const bool /*omit_zeroing_entries*/)
+    Vector<Number, MemorySpace>::reinit(const IndexSet &parallel_partitioner,
+                                        const MPI_Comm  communicator,
+                                        const bool /*omit_zeroing_entries*/)
     {
       vector.reset();
       nonlocal_vector.reset();
@@ -152,12 +152,12 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::reinit(const IndexSet &locally_owned_entries,
-                           const IndexSet &ghost_entries,
-                           const MPI_Comm  communicator,
-                           const bool      vector_writable)
+    Vector<Number, MemorySpace>::reinit(const IndexSet &locally_owned_entries,
+                                        const IndexSet &ghost_entries,
+                                        const MPI_Comm  communicator,
+                                        const bool      vector_writable)
     {
       // release memory before reallocation
       nonlocal_vector.reset();
@@ -201,10 +201,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::reinit(const Vector<Number> &V,
-                           const bool            omit_zeroing_entries)
+    Vector<Number, MemorySpace>::reinit(const Vector<Number, MemorySpace> &V,
+                                        const bool omit_zeroing_entries)
     {
       reinit(V.locally_owned_elements(),
              V.get_mpi_communicator(),
@@ -213,9 +213,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::extract_subvector_to(
+    Vector<Number, MemorySpace>::extract_subvector_to(
       const ArrayView<const types::global_dof_index> &indices,
       ArrayView<Number>                              &elements) const
     {
@@ -271,9 +271,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator=(const Vector<Number> &V)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator=(const Vector<Number, MemorySpace> &V)
     {
       // Distinguish three cases:
       //  - First case: both vectors have the same layout.
@@ -324,9 +324,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator=(const Number s)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator=(const Number s)
     {
       (void)s;
       Assert(s == Number(0.0),
@@ -346,9 +346,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::import_elements(
+    Vector<Number, MemorySpace>::import_elements(
       const ReadWriteVector<Number> &V,
       VectorOperation::values        operation,
       const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
@@ -423,9 +423,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::import_elements(
+    Vector<Number, MemorySpace>::import_elements(
       const ReadWriteVector<Number> &V,
       VectorOperation::values        operation,
       const std::shared_ptr<const Utilities::MPI::CommunicationPatternBase> &)
@@ -435,10 +435,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::import_elements(const ReadWriteVector<Number> &V,
-                                    VectorOperation::values        operation)
+    Vector<Number, MemorySpace>::import_elements(
+      const ReadWriteVector<Number> &V,
+      VectorOperation::values        operation)
     {
       // Create an empty CommunicationPattern
       const Teuchos::RCP<const Utilities::MPI::CommunicationPatternBase>
@@ -449,9 +450,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator*=(const Number factor)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator*=(const Number factor)
     {
       AssertIsFinite(factor);
       vector->scale(factor);
@@ -461,9 +462,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator/=(const Number factor)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator/=(const Number factor)
     {
       AssertIsFinite(factor);
       Assert(factor != Number(0.), ExcZero());
@@ -474,9 +475,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator+=(const Vector<Number> &V)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator+=(
+      const Vector<Number, MemorySpace> &V)
     {
       // If the maps are the same we can update right away.
       if (vector->getMap()->isSameAs(*(V.trilinos_vector().getMap())))
@@ -504,9 +506,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    Vector<Number> &
-    Vector<Number>::operator-=(const Vector<Number> &V)
+    template <typename Number, typename MemorySpace>
+    Vector<Number, MemorySpace> &
+    Vector<Number, MemorySpace>::operator-=(
+      const Vector<Number, MemorySpace> &V)
     {
       this->add(-1., V);
 
@@ -515,9 +518,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Number
-    Vector<Number>::operator*(const Vector<Number> &V) const
+    Vector<Number, MemorySpace>::operator*(
+      const Vector<Number, MemorySpace> &V) const
     {
       Assert(this->size() == V.size(),
              ExcDimensionMismatch(this->size(), V.size()));
@@ -530,9 +534,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Number
-    Vector<Number>::operator()(const size_type index) const
+    Vector<Number, MemorySpace>::operator()(const size_type index) const
     {
       // Get the local index
       const TrilinosWrappers::types::int_type local_index =
@@ -569,9 +573,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::add(const Number a)
+    Vector<Number, MemorySpace>::add(const Number a)
     {
       AssertIsFinite(a);
 
@@ -604,9 +608,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::add(const Number a, const Vector<Number> &V)
+    Vector<Number, MemorySpace>::add(const Number                       a,
+                                     const Vector<Number, MemorySpace> &V)
     {
       AssertIsFinite(a);
 
@@ -622,12 +627,12 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::add(const Number          a,
-                        const Vector<Number> &V,
-                        const Number          b,
-                        const Vector<Number> &W)
+    Vector<Number, MemorySpace>::add(const Number                       a,
+                                     const Vector<Number, MemorySpace> &V,
+                                     const Number                       b,
+                                     const Vector<Number, MemorySpace> &W)
     {
       AssertIsFinite(a);
       AssertIsFinite(b);
@@ -646,11 +651,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::sadd(const Number          s,
-                         const Number          a,
-                         const Vector<Number> &V)
+    Vector<Number, MemorySpace>::sadd(const Number                       s,
+                                      const Number                       a,
+                                      const Vector<Number, MemorySpace> &V)
     {
       AssertIsFinite(s);
       AssertIsFinite(a);
@@ -661,16 +666,17 @@ namespace LinearAlgebra
 
       *this *= s;
 
-      Vector<Number> tmp(V);
+      Vector<Number, MemorySpace> tmp(V);
       tmp *= a;
       *this += tmp;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::scale(const Vector<Number> &scaling_factors)
+    Vector<Number, MemorySpace>::scale(
+      const Vector<Number, MemorySpace> &scaling_factors)
     {
       Assert(vector->getMap()->isSameAs(
                *(scaling_factors.trilinos_vector().getMap())),
@@ -685,9 +691,10 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::equ(const Number a, const Vector<Number> &V)
+    Vector<Number, MemorySpace>::equ(const Number                       a,
+                                     const Vector<Number, MemorySpace> &V)
     {
       AssertIsFinite(a);
 
@@ -707,9 +714,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     bool
-    Vector<Number>::all_zero() const
+    Vector<Number, MemorySpace>::all_zero() const
     {
       // get a representation of the vector and
       // loop over all the elements
@@ -738,9 +745,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Number
-    Vector<Number>::mean_value() const
+    Vector<Number, MemorySpace>::mean_value() const
     {
       Assert(!has_ghost_elements(), ExcGhostsPresent());
 
@@ -749,9 +756,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    typename Vector<Number>::real_type
-    Vector<Number>::l1_norm() const
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::real_type
+    Vector<Number, MemorySpace>::l1_norm() const
     {
       Assert(!has_ghost_elements(), ExcGhostsPresent());
 
@@ -760,9 +767,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    typename Vector<Number>::real_type
-    Vector<Number>::l2_norm() const
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::real_type
+    Vector<Number, MemorySpace>::l2_norm() const
     {
       Assert(!has_ghost_elements(), ExcGhostsPresent());
 
@@ -771,9 +778,9 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    typename Vector<Number>::real_type
-    Vector<Number>::linfty_norm() const
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::real_type
+    Vector<Number, MemorySpace>::linfty_norm() const
     {
       Assert(!has_ghost_elements(), ExcGhostsPresent());
 
@@ -782,11 +789,12 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Number
-    Vector<Number>::add_and_dot(const Number          a,
-                                const Vector<Number> &V,
-                                const Vector<Number> &W)
+    Vector<Number, MemorySpace>::add_and_dot(
+      const Number                       a,
+      const Vector<Number, MemorySpace> &V,
+      const Vector<Number, MemorySpace> &W)
     {
       AssertIsFinite(a);
 
@@ -797,27 +805,27 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
-    typename Vector<Number>::size_type
-    Vector<Number>::size() const
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::size_type
+    Vector<Number, MemorySpace>::size() const
     {
       return vector->getGlobalLength();
     }
 
 
 
-    template <typename Number>
-    typename Vector<Number>::size_type
-    Vector<Number>::locally_owned_size() const
+    template <typename Number, typename MemorySpace>
+    typename Vector<Number, MemorySpace>::size_type
+    Vector<Number, MemorySpace>::locally_owned_size() const
     {
       return vector->getLocalLength();
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     MPI_Comm
-    Vector<Number>::get_mpi_communicator() const
+    Vector<Number, MemorySpace>::get_mpi_communicator() const
     {
       return Utilities::Trilinos::teuchos_comm_to_mpi_comm(
         vector->getMap()->getComm());
@@ -825,18 +833,19 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     IndexSet
-    Vector<Number>::locally_owned_elements() const
+    Vector<Number, MemorySpace>::locally_owned_elements() const
     {
       return IndexSet(vector->getMap());
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::compress(const VectorOperation::values operation)
+    Vector<Number, MemorySpace>::compress(
+      const VectorOperation::values operation)
     {
       Assert(has_ghost == false,
              ExcMessage(
@@ -866,49 +875,49 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     const Tpetra::Vector<Number, int, types::signed_global_dof_index> &
-    Vector<Number>::trilinos_vector() const
+    Vector<Number, MemorySpace>::trilinos_vector() const
     {
       return *vector;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Tpetra::Vector<Number, int, types::signed_global_dof_index> &
-    Vector<Number>::trilinos_vector()
+    Vector<Number, MemorySpace>::trilinos_vector()
     {
       return *vector;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Teuchos::RCP<Tpetra::Vector<Number, int, types::signed_global_dof_index>>
-    Vector<Number>::trilinos_rcp()
+    Vector<Number, MemorySpace>::trilinos_rcp()
     {
       return vector;
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     Teuchos::RCP<
       const Tpetra::Vector<Number, int, types::signed_global_dof_index>>
-    Vector<Number>::trilinos_rcp() const
+    Vector<Number, MemorySpace>::trilinos_rcp() const
     {
       return vector.getConst();
     }
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::print(std::ostream      &out,
-                          const unsigned int precision,
-                          const bool         scientific,
-                          const bool         across) const
+    Vector<Number, MemorySpace>::print(std::ostream      &out,
+                                       const unsigned int precision,
+                                       const bool         scientific,
+                                       const bool         across) const
     {
       AssertThrow(out.fail() == false, ExcIO());
       boost::io::ios_flags_saver restore_flags(out);
@@ -961,18 +970,18 @@ namespace LinearAlgebra
     }
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     MPI_Comm
-    Vector<Number>::mpi_comm() const
+    Vector<Number, MemorySpace>::mpi_comm() const
     {
       return Utilities::Trilinos::teuchos_comm_to_mpi_comm(
         vector->getMap()->getComm());
     }
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     std::size_t
-    Vector<Number>::memory_consumption() const
+    Vector<Number, MemorySpace>::memory_consumption() const
     {
       return sizeof(*this) +
              vector->getLocalLength() *
@@ -981,10 +990,11 @@ namespace LinearAlgebra
 
 
 
-    template <typename Number>
+    template <typename Number, typename MemorySpace>
     void
-    Vector<Number>::create_tpetra_comm_pattern(const IndexSet &source_index_set,
-                                               const MPI_Comm  mpi_comm)
+    Vector<Number, MemorySpace>::create_tpetra_comm_pattern(
+      const IndexSet &source_index_set,
+      const MPI_Comm  mpi_comm)
     {
       source_stored_elements = source_index_set;