From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Mon, 19 Feb 2024 23:10:50 +0000 (-0700)
Subject: Add Tpetra SparseMatrix::matrix_scalar_product() and matrix_norm_square().
X-Git-Tag: relicensing~17^2~2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=a0b947e0e4fd0caeb3d60cb035213503322bee5b;p=dealii.git

Add Tpetra SparseMatrix::matrix_scalar_product() and matrix_norm_square().
---

diff --git a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
index ad91ed457a..4ca88bde7e 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.h
@@ -784,6 +784,53 @@ namespace LinearAlgebra
       Tvmult_add(Vector<Number, MemorySpace>       &dst,
                  const Vector<Number, MemorySpace> &src) const;
 
+      /**
+       * Return the square of the norm of the vector $v$ with respect to the
+       * norm induced by this matrix, i.e., $\left(v,Mv\right)$. This is useful,
+       * e.g. in the finite element context, where the $L_2$ norm of a function
+       * equals the matrix norm with respect to the @ref GlossMassMatrix "mass matrix" of the vector
+       * representing the nodal values of the finite element function.
+       *
+       * Obviously, the matrix needs to be quadratic for this operation.
+       *
+       * The implementation of this function is not as efficient as the one in
+       * the @p SparseMatrix class used in deal.II (i.e. the original one, not
+       * the Trilinos wrapper class) since Trilinos doesn't support this
+       * operation and needs a temporary vector.
+       *
+       * The vector has to be initialized with the same IndexSet the matrix
+       * was initialized with.
+       *
+       * In case of a localized Vector, this function will only work when
+       * running on one processor, since the matrix object is inherently
+       * distributed. Otherwise, an exception will be thrown.
+       */
+      Number
+      matrix_norm_square(const Vector<Number, MemorySpace> &v) const;
+
+      /**
+       * Compute the matrix scalar product $\left(u,Mv\right)$.
+       *
+       * The implementation of this function is not as efficient as the one in
+       * the @p SparseMatrix class used in deal.II (i.e. the original one, not
+       * the Trilinos wrapper class) since Trilinos doesn't support this
+       * operation and needs a temporary vector.
+       *
+       * The vector @p u has to be initialized with the same IndexSet that
+       * was used for the row indices of the matrix and the vector @p v has
+       * to be initialized with the same IndexSet that was used for the
+       * column indices of the matrix.
+       *
+       * In case of a localized Vector, this function will only work when
+       * running on one processor, since the matrix object is inherently
+       * distributed. Otherwise, an exception will be thrown.
+       *
+       * This function is only implemented for square matrices.
+       */
+      Number
+      matrix_scalar_product(const Vector<Number, MemorySpace> &u,
+                            const Vector<Number, MemorySpace> &v) const;
+
       /**
        * Compute the residual of an equation <i>Mx=b</i>, where the residual is
        * defined to be <i>r=b-Mx</i>. Write the residual into @p dst. 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 a8c1497817..9c20d52713 100644
--- a/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
+++ b/include/deal.II/lac/trilinos_tpetra_sparse_matrix.templates.h
@@ -801,6 +801,43 @@ namespace LinearAlgebra
 
 
 
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::matrix_norm_square(
+      const Vector<Number, MemorySpace> &v) const
+    {
+      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
+      Assert(matrix->getRowMap()->isSameAs(*matrix->getDomainMap()),
+             ExcNotQuadratic());
+
+      Vector<Number, MemorySpace> temp_vector;
+      temp_vector.reinit(v, true);
+
+      vmult(temp_vector, v);
+      return temp_vector * v;
+    }
+
+
+
+    template <typename Number, typename MemorySpace>
+    Number
+    SparseMatrix<Number, MemorySpace>::matrix_scalar_product(
+      const Vector<Number, MemorySpace> &u,
+      const Vector<Number, MemorySpace> &v) const
+    {
+      Assert(matrix->isFillComplete(), ExcMatrixNotCompressed());
+      Assert(matrix->getRowMap()->isSameAs(*matrix->getDomainMap()),
+             ExcNotQuadratic());
+
+      Vector<Number, MemorySpace> temp_vector;
+      temp_vector.reinit(v, true);
+
+      vmult(temp_vector, v);
+      return u * temp_vector;
+    }
+
+
+
     template <typename Number, typename MemorySpace>
     Number
     SparseMatrix<Number, MemorySpace>::frobenius_norm() const