From: Wolfgang Bangerth 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 &dst, const Vector &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 &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 &u, + const Vector &v) const; + /** * Compute the residual of an equation Mx=b, where the residual is * defined to be r=b-Mx. 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 + Number + SparseMatrix::matrix_norm_square( + const Vector &v) const + { + Assert(matrix->isFillComplete(), ExcMatrixNotCompressed()); + Assert(matrix->getRowMap()->isSameAs(*matrix->getDomainMap()), + ExcNotQuadratic()); + + Vector temp_vector; + temp_vector.reinit(v, true); + + vmult(temp_vector, v); + return temp_vector * v; + } + + + + template + Number + SparseMatrix::matrix_scalar_product( + const Vector &u, + const Vector &v) const + { + Assert(matrix->isFillComplete(), ExcMatrixNotCompressed()); + Assert(matrix->getRowMap()->isSameAs(*matrix->getDomainMap()), + ExcNotQuadratic()); + + Vector temp_vector; + temp_vector.reinit(v, true); + + vmult(temp_vector, v); + return u * temp_vector; + } + + + template Number SparseMatrix::frobenius_norm() const