From c94cb9d961bd86969b5912cb691eb3fdecb5fe29 Mon Sep 17 00:00:00 2001 From: Jie Cheng Date: Wed, 20 Dec 2017 15:47:56 -0500 Subject: [PATCH] Add PETScWrappers::MatrixBase::mmult --- include/deal.II/lac/petsc_matrix_base.h | 47 ++++++++++++- source/lac/petsc_matrix_base.cc | 88 +++++++++++++++++++++++++ 2 files changed, 133 insertions(+), 2 deletions(-) diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h index 72cb897e30..7087cdeeb4 100644 --- a/include/deal.II/lac/petsc_matrix_base.h +++ b/include/deal.II/lac/petsc_matrix_base.h @@ -26,6 +26,7 @@ # include # include # include +# include # include @@ -41,7 +42,6 @@ template class BlockMatrixBase; namespace PETScWrappers { // forward declarations - class VectorBase; class MatrixBase; namespace MatrixIterators @@ -784,7 +784,6 @@ namespace PETScWrappers void Tvmult_add (VectorBase &dst, const VectorBase &src) 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 @@ -840,6 +839,13 @@ namespace PETScWrappers */ operator Mat () const; + /** + * Return a reference to the underlying PETSc type. It can be used to + * modify the underlying data, so use it only when you know what you + * are doing. + */ + Mat &petsc_matrix (); + /** * Make an in-place transpose of a matrix. */ @@ -952,7 +958,44 @@ namespace PETScWrappers */ void prepare_set(); + /** + * Base function to perform the matrix-matrix multiplication $C = AB$, + * or, if a vector $V$ whose size is compatible with B is given, + * $C = A \text{diag}(V) B$, where $\text{diag}(V)$ defines a + * diagonal matrix with the vector entries. + * + * This function assumes that the calling matrix $A$ and $B$ + * have compatible sizes. The size of $C$ will be set within this + * function. + * + * The content as well as the sparsity pattern of the matrix $C$ will be + * reset by this function, so make sure that the sparsity pattern is not + * used somewhere else in your program. This is an expensive operation, so + * think twice before you use this function. + */ + void mmult (MatrixBase &C, + const MatrixBase &B, + const VectorBase &V = VectorBase()) const; + /** + * Base function to perform the matrix-matrix multiplication with + * the transpose of this, i.e., $C = A^T B$, or, + * if an optional vector $V$ whose size is compatible with $B$ is given, + * $C = A^T \text{diag}(V) B$, where $\text{diag}(V)$ defines a + * diagonal matrix with the vector entries. + * + * This function assumes that the calling matrix $A$ and $B$ + * have compatible sizes. The size of $C$ will be set within this + * function. + * + * The content as well as the sparsity pattern of the matrix $C$ will be + * changed by this function, so make sure that the sparsity pattern is not + * used somewhere else in your program. This is an expensive operation, so + * think twice before you use this function. + */ + void Tmmult (MatrixBase &C, + const MatrixBase &B, + const VectorBase &V = VectorBase()) const; private: diff --git a/source/lac/petsc_matrix_base.cc b/source/lac/petsc_matrix_base.cc index efc57ad532..cd1e821f99 100644 --- a/source/lac/petsc_matrix_base.cc +++ b/source/lac/petsc_matrix_base.cc @@ -489,6 +489,89 @@ namespace PETScWrappers } + namespace internals + { + void perform_mmult (const MatrixBase &inputleft, + const MatrixBase &inputright, + MatrixBase &result, + const VectorBase &V, + const bool transpose_left) + { + const bool use_vector = (V.size() == inputright.m() ? true : false); + if (transpose_left == false) + { + Assert (inputleft.n() == inputright.m(), + ExcDimensionMismatch(inputleft.n(), inputright.m())); + } + else + { + Assert (inputleft.m() == inputright.m(), + ExcDimensionMismatch(inputleft.m(), inputright.m())); + } + + result.clear(); + + PetscErrorCode ierr; + + if (use_vector == false) + { + if (transpose_left) + { + ierr = MatTransposeMatMult (inputleft, + inputright, + MAT_INITIAL_MATRIX, + PETSC_DEFAULT, + &result.petsc_matrix()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + else + { + ierr = MatMatMult (inputleft, + inputright, + MAT_INITIAL_MATRIX, + PETSC_DEFAULT, + &result.petsc_matrix()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + } + else + { + Mat tmp; + ierr = MatDuplicate (inputleft, MAT_COPY_VALUES, &tmp); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + if (transpose_left) + { + ierr = MatTranspose(tmp, MAT_REUSE_MATRIX, &tmp); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + ierr = MatDiagonalScale (tmp, NULL, V); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + ierr = MatMatMult (tmp, + inputright, + MAT_INITIAL_MATRIX, + PETSC_DEFAULT, + &result.petsc_matrix()); + AssertThrow (ierr == 0, ExcPETScError(ierr)); + } + } + } + + void + MatrixBase::mmult (MatrixBase &C, + const MatrixBase &B, + const VectorBase &V) const + { + internals::perform_mmult (*this, B, C, V, false); + } + + void + MatrixBase::Tmmult (MatrixBase &C, + const MatrixBase &B, + const VectorBase &V) const + { + internals::perform_mmult (*this, B, C, V, true); + } + PetscScalar MatrixBase::residual (VectorBase &dst, const VectorBase &x, @@ -511,6 +594,11 @@ namespace PETScWrappers return matrix; } + Mat &MatrixBase::petsc_matrix () + { + return matrix; + } + void MatrixBase::transpose () { -- 2.39.5