From: Sebastian Kinnewig <kinnewig@ifam.uni-hannover.de>
Date: Fri, 9 Feb 2024 10:40:23 +0000 (+0100)
Subject: Template TpetraWrappers::SparseMatrix::vmult on the vector type.
X-Git-Tag: v9.6.0-rc1~518^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=9db41cde025579225e9c9d2f730d6f9d3f672217;p=dealii.git

Template TpetraWrappers::SparseMatrix::vmult on the vector type.
---

diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
index f29536400e..e466f89851 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
@@ -27,11 +27,14 @@
 #  include <deal.II/lac/sparsity_pattern.h>
 #  include <deal.II/lac/trilinos_tpetra_sparsity_pattern.h>
 #  include <deal.II/lac/trilinos_tpetra_vector.h>
+#  include <deal.II/lac/vector.h>
 
 // Tpetra includes
 #  include <Tpetra_Core.hpp>
 #  include <Tpetra_CrsMatrix.hpp>
 
+#  include <type_traits>
+
 
 DEAL_II_NAMESPACE_OPEN
 
@@ -162,6 +165,12 @@ namespace LinearAlgebra
       using GraphType =
         Tpetra::CrsGraph<int, dealii::types::signed_global_dof_index, NodeType>;
 
+      /**
+       * Typedef for Tpetra::Vector
+       */
+      using VectorType = Tpetra::
+        Vector<Number, int, dealii::types::signed_global_dof_index, NodeType>;
+
       /**
        * @name Constructors and initialization.
        */
@@ -799,9 +808,9 @@ namespace LinearAlgebra
        * initialized with the same IndexSet that was used for the column indices
        * of the matrix.
        */
+      template <typename InputVectorType>
       void
-      vmult(Vector<Number, MemorySpace>       &dst,
-            const Vector<Number, MemorySpace> &src) const;
+      vmult(InputVectorType &dst, const InputVectorType &src) const;
 
       /*
        * Matrix-vector multiplication: let <i>dst = M<sup>T</sup>*src</i> with
@@ -810,9 +819,9 @@ namespace LinearAlgebra
        *
        * Source and destination must not be the same vector.
        */
+      template <typename InputVectorType>
       void
-      Tvmult(Vector<Number, MemorySpace>       &dst,
-             const Vector<Number, MemorySpace> &src) const;
+      Tvmult(InputVectorType &dst, const InputVectorType &src) const;
 
       /**
        * Adding matrix-vector multiplication. Add <i>M*src</i> on <i>dst</i>
@@ -820,10 +829,9 @@ namespace LinearAlgebra
        *
        * Source and destination must not be the same vector.
        */
+      template <typename InputVectorType>
       void
-      vmult_add(Vector<Number, MemorySpace>       &dst,
-                const Vector<Number, MemorySpace> &src) const;
-
+      vmult_add(InputVectorType &dst, const InputVectorType &src) const;
 
       /**
        * Adding matrix-vector multiplication. Add <i>M<sup>T</sup>*src</i> to
@@ -832,9 +840,9 @@ namespace LinearAlgebra
        *
        * Source and destination must not be the same vector.
        */
+      template <typename InputVectorType>
       void
-      Tvmult_add(Vector<Number, MemorySpace>       &dst,
-                 const Vector<Number, MemorySpace> &src) const;
+      Tvmult_add(InputVectorType &dst, const InputVectorType &src) const;
 
       /**
        * Return the square of the norm of the vector $v$ with respect to the
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 e0660c7729..4a49d36982 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -31,6 +31,114 @@ namespace LinearAlgebra
 
   namespace TpetraWrappers
   {
+    namespace internal
+    {
+      template <typename Number, typename MemorySpace>
+      void
+      apply(const SparseMatrix<Number, MemorySpace> &M,
+            const Vector<Number, MemorySpace>       &src,
+            Vector<Number, MemorySpace>             &dst,
+            Teuchos::ETransp                         mode = Teuchos::NO_TRANS,
+            Number alpha = Teuchos::ScalarTraits<Number>::one(),
+            Number beta  = Teuchos::ScalarTraits<Number>::zero())
+      {
+        Assert(&src != &dst,
+               SparseMatrix<double>::ExcSourceEqualsDestination());
+        Assert(M.trilinos_matrix().isFillComplete(),
+               SparseMatrix<double>::ExcMatrixNotCompressed());
+
+        if (mode == Teuchos::NO_TRANS)
+          {
+            Assert(src.trilinos_vector().getMap()->isSameAs(
+                     *M.trilinos_matrix().getDomainMap()),
+                   SparseMatrix<double>::ExcColMapMissmatch());
+            Assert(dst.trilinos_vector().getMap()->isSameAs(
+                     *M.trilinos_matrix().getRangeMap()),
+                   SparseMatrix<double>::ExcDomainMapMissmatch());
+          }
+        else
+          {
+            Assert(dst.trilinos_vector().getMap()->isSameAs(
+                     *M.trilinos_matrix().getDomainMap()),
+                   SparseMatrix<double>::ExcColMapMissmatch());
+            Assert(src.trilinos_vector().getMap()->isSameAs(
+                     *M.trilinos_matrix().getRangeMap()),
+                   SparseMatrix<double>::ExcDomainMapMissmatch());
+          }
+
+        M.trilinos_matrix().apply(
+          src.trilinos_vector(), dst.trilinos_vector(), mode, alpha, beta);
+      }
+
+
+
+      template <typename Number, typename MemorySpace>
+      void
+      apply(const SparseMatrix<Number, MemorySpace> &M,
+            const dealii::Vector<Number>            &src,
+            dealii::Vector<Number>                  &dst,
+            Teuchos::ETransp                         mode = Teuchos::NO_TRANS,
+            Number alpha = Teuchos::ScalarTraits<Number>::one(),
+            Number beta  = Teuchos::ScalarTraits<Number>::zero())
+      {
+        Assert(&src != &dst,
+               SparseMatrix<double>::ExcSourceEqualsDestination());
+        Assert(M.trilinos_matrix().isFillComplete(),
+               SparseMatrix<double>::ExcMatrixNotCompressed());
+
+        // get the size of the input vectors:
+        const size_type dst_local_size = dst.end() - dst.begin();
+        const size_type src_local_size = src.end() - src.begin();
+
+        // For the dst vector:
+        Kokkos::View<Number **, Kokkos::LayoutLeft, Kokkos::HostSpace>
+          kokkos_view_dst(dst.begin(), dst_local_size, 1);
+
+        // get a Kokkos::DualView
+        auto mirror_view_dst = Kokkos::create_mirror_view_and_copy(
+          typename MemorySpace::kokkos_space{}, kokkos_view_dst);
+        typename SparseMatrix<Number, MemorySpace>::VectorType::dual_view_type
+          kokkos_dual_view_dst(mirror_view_dst, kokkos_view_dst);
+
+        // create the Tpetra::Vector
+        typename SparseMatrix<Number, MemorySpace>::VectorType tpetra_dst(
+          M.trilinos_matrix().getRangeMap(), kokkos_dual_view_dst);
+
+        // For the src vector:
+        // create a Kokkos::View from the src vector
+        Kokkos::View<Number **, Kokkos::LayoutLeft, Kokkos::HostSpace>
+          kokkos_view_src(const_cast<Number *>(src.begin()), src_local_size, 1);
+
+        // get a Kokkos::DualView
+        auto mirror_view_src = Kokkos::create_mirror_view_and_copy(
+          typename MemorySpace::kokkos_space{}, kokkos_view_src);
+        typename SparseMatrix<Number, MemorySpace>::VectorType::dual_view_type
+          kokkos_dual_view_src(mirror_view_src, kokkos_view_src);
+
+        // create the Tpetra::Vector
+        typename SparseMatrix<Number, MemorySpace>::VectorType tpetra_src(
+          M.trilinos_matrix().getDomainMap(), kokkos_dual_view_src);
+
+        M.trilinos_matrix().apply(tpetra_src, tpetra_dst, mode, alpha, beta);
+      }
+
+
+
+      template <typename Number, typename MemorySpace, typename VectorType>
+      void
+      apply(SparseMatrix<Number, MemorySpace> &,
+            const VectorType &,
+            VectorType &,
+            Teuchos::ETransp,
+            Number,
+            Number)
+      {
+        DEAL_II_NOT_IMPLEMENTED();
+      }
+    } // namespace internal
+
+
+
     // reinit_matrix():
     namespace
     {
@@ -1170,81 +1278,58 @@ namespace LinearAlgebra
 
 
     // Multiplications
-
     template <typename Number, typename MemorySpace>
+    template <typename InputVectorType>
     void
-    SparseMatrix<Number, MemorySpace>::vmult(
-      Vector<Number, MemorySpace>       &dst,
-      const Vector<Number, MemorySpace> &src) const
+    SparseMatrix<Number, MemorySpace>::vmult(InputVectorType       &dst,
+                                             const InputVectorType &src) const
     {
-      Assert(&src != &dst, ExcSourceEqualsDestination());
-      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
-      Assert(src.trilinos_vector().getMap()->isSameAs(*matrix->getDomainMap()),
-             ExcColMapMissmatch());
-      Assert(dst.trilinos_vector().getMap()->isSameAs(*matrix->getRangeMap()),
-             ExcDomainMapMissmatch());
-      matrix->apply(src.trilinos_vector(), dst.trilinos_vector());
+      internal::apply(*this, src, dst);
     }
 
 
 
     template <typename Number, typename MemorySpace>
+    template <typename InputVectorType>
     void
-    SparseMatrix<Number, MemorySpace>::Tvmult(
-      Vector<Number, MemorySpace>       &dst,
-      const Vector<Number, MemorySpace> &src) const
+    SparseMatrix<Number, MemorySpace>::Tvmult(InputVectorType       &dst,
+                                              const InputVectorType &src) const
     {
-      Assert(&src != &dst, ExcSourceEqualsDestination());
-      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
-      Assert(dst.trilinos_vector().getMap()->isSameAs(*matrix->getDomainMap()),
-             ExcColMapMissmatch());
-      Assert(src.trilinos_vector().getMap()->isSameAs(*matrix->getRangeMap()),
-             ExcDomainMapMissmatch());
-      matrix->apply(src.trilinos_vector(),
-                    dst.trilinos_vector(),
-                    Teuchos::TRANS);
+      internal::apply(*this, src, dst, Teuchos::TRANS);
     }
 
 
 
     template <typename Number, typename MemorySpace>
+    template <typename InputVectorType>
     void
     SparseMatrix<Number, MemorySpace>::vmult_add(
-      Vector<Number, MemorySpace>       &dst,
-      const Vector<Number, MemorySpace> &src) const
+      InputVectorType       &dst,
+      const InputVectorType &src) const
     {
-      Assert(&src != &dst, ExcSourceEqualsDestination());
-      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
-      Assert(src.trilinos_vector().getMap()->isSameAs(*matrix->getDomainMap()),
-             ExcColMapMissmatch());
-      Assert(dst.trilinos_vector().getMap()->isSameAs(*matrix->getRangeMap()),
-             ExcDomainMapMissmatch());
-      matrix->apply(src.trilinos_vector(),
-                    dst.trilinos_vector(),
-                    Teuchos::NO_TRANS,
-                    Teuchos::ScalarTraits<Number>::one(),
-                    Teuchos::ScalarTraits<Number>::one());
+      internal::apply(*this,
+                      src,
+                      dst,
+                      Teuchos::NO_TRANS,
+                      Teuchos::ScalarTraits<Number>::one(),
+                      Teuchos::ScalarTraits<Number>::one());
     }
 
 
 
     template <typename Number, typename MemorySpace>
+    template <typename InputVectorType>
     void
     SparseMatrix<Number, MemorySpace>::Tvmult_add(
-      Vector<Number, MemorySpace>       &dst,
-      const Vector<Number, MemorySpace> &src) const
+      InputVectorType       &dst,
+      const InputVectorType &src) const
     {
-      Assert(&src != &dst, ExcSourceEqualsDestination());
-      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
-      Assert(dst.trilinos_vector().getMap()->isSameAs(*matrix->getDomainMap()),
-             ExcColMapMissmatch());
-      Assert(src.trilinos_vector().getMap()->isSameAs(*matrix->getRangeMap()),
-             ExcDomainMapMissmatch());
-      matrix->apply(src.trilinos_vector(),
-                    dst.trilinos_vector(),
-                    Teuchos::TRANS,
-                    Teuchos::ScalarTraits<Number>::one(),
-                    Teuchos::ScalarTraits<Number>::one());
+      internal::apply(*this,
+                      src,
+                      dst,
+                      Teuchos::TRANS,
+                      Teuchos::ScalarTraits<Number>::one(),
+                      Teuchos::ScalarTraits<Number>::one());
     }
 
 
diff --git a/source/lac/trilinos_tpetra_sparse_matrix.cc b/source/lac/trilinos_tpetra_sparse_matrix.cc
index 76c5335856..3cc43d6226 100644
--- a/source/lac/trilinos_tpetra_sparse_matrix.cc
+++ b/source/lac/trilinos_tpetra_sparse_matrix.cc
@@ -46,6 +46,37 @@ namespace LinearAlgebra
     template void
     SparseMatrix<double>::reinit(const dealii::DynamicSparsityPattern &);
 
+    template void
+    SparseMatrix<double>::vmult(Vector<double>       &dst,
+                                const Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::Tvmult(Vector<double>       &dst,
+                                 const Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::vmult_add(Vector<double>       &dst,
+                                    const Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::Tvmult_add(Vector<double>       &dst,
+                                     const Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::vmult(::dealii::Vector<double>       &dst,
+                                const ::dealii::Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::Tvmult(::dealii::Vector<double>       &dst,
+                                 const ::dealii::Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::vmult_add(::dealii::Vector<double>       &dst,
+                                    const ::dealii::Vector<double> &src) const;
+
+    template void
+    SparseMatrix<double>::Tvmult_add(::dealii::Vector<double>       &dst,
+                                     const ::dealii::Vector<double> &src) const;
   } // namespace TpetraWrappers
 } // namespace LinearAlgebra
 #  endif // DOXYGEN