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 ad91ed457a65947698d4e6c4b84ec94309900497..4ca88bde7e59e0ad9e1b2fc4a12c52944e6239b3 100644 (file)
--- 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 a8c14978179dc610afc7809031441e84dd40b254..9c20d527134bc58b25f425b71fb8650289f76bcd 100644 (file)
--- 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