Add PETScWrappers::MatrixBase::mmult

diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h
index 72cb897e305f04636d4845f4d444285ce81dbb6c..7087cdeeb4bf473abc2f8b9fe0d3625e929a0946 100644 (file)
--- a/include/deal.II/lac/petsc_matrix_base.h
+++ b/include/deal.II/lac/petsc_matrix_base.h
@@ -26,6 +26,7 @@
 #  include <deal.II/lac/full_matrix.h>
 #  include <deal.II/lac/petsc_compatibility.h>
 #  include <deal.II/lac/vector_operation.h>
+#  include <deal.II/lac/petsc_vector_base.h>
 
 #  include <petscmat.h>
 
@@ -41,7 +42,6 @@ template <typename Matrix> 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 <i>Mx=b</i>, where the residual is
      * defined to be <i>r=b-Mx</i>. 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 <tt>this</tt>, 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 efc57ad532c53d44020eda34944ffd2239d567d9..cd1e821f99722e8f4285a21aba87b6d66259b0e2 100644 (file)
--- 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 ()
   {