PETScWrappers:BlockSparseMatrix: constructor from PETSc Mat

diff --git a/include/deal.II/lac/petsc_block_sparse_matrix.h b/include/deal.II/lac/petsc_block_sparse_matrix.h
index f62d718b2ab17fbfc64f05be287785c1dc188c56..626974e70a0e2775ecf722dae72ba61bcabae3ad 100644 (file)
--- a/include/deal.II/lac/petsc_block_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_block_sparse_matrix.h
@@ -102,6 +102,11 @@ namespace PETScWrappers
        */
       BlockSparseMatrix() = default;
 
+      /**
+       * Create a BlockSparseMatrix with a PETSc Mat
+       */
+      explicit BlockSparseMatrix(const Mat &);
+
       /**
        * Destructor.
        */
@@ -289,6 +294,12 @@ namespace PETScWrappers
       Mat &
       petsc_matrix();
 
+      /**
+       * This method assigns the PETSc Mat to the instance of the class.
+       */
+      void
+      assign_petsc_matrix(Mat A);
+
     private:
       Mat petsc_nest_matrix = nullptr;
     };
@@ -299,6 +310,12 @@ namespace PETScWrappers
 
     // ------------- inline and template functions -----------------
 
+    inline BlockSparseMatrix::BlockSparseMatrix(const Mat &A)
+      : BlockSparseMatrix()
+    {
+      this->assign_petsc_matrix(A);
+    }
+
     inline BlockSparseMatrix &
     BlockSparseMatrix::operator=(const double d)
     {

diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h
index 2fd3d37f0c1f26318a7bf6fe79f9048a583e0975..366d4172eedaf0ffadd59535cfb557a38c29cdfa 100644 (file)
--- a/include/deal.II/lac/petsc_matrix_base.h
+++ b/include/deal.II/lac/petsc_matrix_base.h
@@ -344,6 +344,13 @@ namespace PETScWrappers
      */
     virtual ~MatrixBase() override;
 
+    /**
+     * This method assignes the PETSc Mat to the instance of the class.
+     *
+     */
+    void
+    assign_petsc_matrix(Mat A);
+
     /**
      * This operator assigns a scalar to a matrix. Since this does usually not
      * make much sense (should we set all matrix entries to this value? Only

diff --git a/include/deal.II/lac/petsc_sparse_matrix.h b/include/deal.II/lac/petsc_sparse_matrix.h
index 08d7734571adcf133f393b75d1095e15800d9d83..a3ced65c7c06bbf24c00d7d434b4f1e37435e9a7 100644 (file)
--- a/include/deal.II/lac/petsc_sparse_matrix.h
+++ b/include/deal.II/lac/petsc_sparse_matrix.h
@@ -75,6 +75,12 @@ namespace PETScWrappers
      */
     SparseMatrix();
 
+    /**
+     * Initialize a Matrix from a PETSc Mat object. Note that we do not copy
+     * the matrix.
+     */
+    explicit SparseMatrix(const Mat &);
+
     /**
      * Create a sparse matrix of dimensions @p m times @p n, with an initial
      * guess of @p n_nonzero_per_row nonzero elements per row. PETSc is able
@@ -397,6 +403,12 @@ namespace PETScWrappers
        */
       SparseMatrix();
 
+      /**
+       * Initialize a SparseMatrix from a PETSc Mat object. Note that we do not
+       * copy the matrix.
+       */
+      explicit SparseMatrix(const Mat &);
+
       /**
        * Destructor to free the PETSc object.
        */

diff --git a/source/lac/petsc_matrix_base.cc b/source/lac/petsc_matrix_base.cc
index b7dd5c688382dfb58c88b3b94b930158eb6372ac..c7363ba9c7aec9afb6a16064392836491e93ba02 100644 (file)
--- a/source/lac/petsc_matrix_base.cc
+++ b/source/lac/petsc_matrix_base.cc
@@ -89,14 +89,23 @@ namespace PETScWrappers
     AssertThrow(ierr == 0, ExcPETScError(ierr));
   }
 
+  void
+  MatrixBase::assign_petsc_matrix(Mat A)
+  {
+    AssertThrow(last_action == ::dealii::VectorOperation::unknown,
+                ExcMessage("Cannot assign a new Mat"));
+    PetscErrorCode ierr =
+      PetscObjectReference(reinterpret_cast<PetscObject>(A));
+    AssertThrow(ierr == 0, ExcPETScError(ierr));
+    destroy_matrix(matrix);
+    matrix = A;
+  }
 
   MatrixBase::~MatrixBase()
   {
     destroy_matrix(matrix);
   }
 
-
-
   void
   MatrixBase::clear()
   {

diff --git a/source/lac/petsc_parallel_block_sparse_matrix.cc b/source/lac/petsc_parallel_block_sparse_matrix.cc
index a290cfb602bc1cc7b9032b2aba699e769a384440..ae50044645599801f0420c2bf86112cabe7f3fe6 100644 (file)
--- a/source/lac/petsc_parallel_block_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_block_sparse_matrix.cc
@@ -188,6 +188,56 @@ namespace PETScWrappers
       return petsc_nest_matrix;
     }
 
+    void
+    BlockSparseMatrix::assign_petsc_matrix(Mat A)
+    {
+      clear();
+
+      PetscBool isnest;
+      PetscInt  nr = 1, nc = 1;
+
+      PetscErrorCode ierr =
+        PetscObjectTypeCompare(reinterpret_cast<PetscObject>(A),
+                               MATNEST,
+                               &isnest);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
+      std::vector<Mat> mats;
+      if (isnest)
+        {
+          ierr = MatNestGetSize(A, &nr, &nc);
+          AssertThrow(ierr == 0, ExcPETScError(ierr));
+          for (PetscInt i = 0; i < nr; ++i)
+            {
+              for (PetscInt j = 0; j < nc; ++j)
+                {
+                  Mat sA;
+                  ierr = MatNestGetSubMat(A, i, j, &sA);
+                  mats.push_back(sA);
+                }
+            }
+        }
+      else
+        {
+          mats.push_back(A);
+        }
+
+      std::vector<size_type> r_block_sizes(nr, 0);
+      std::vector<size_type> c_block_sizes(nc, 0);
+      this->row_block_indices.reinit(r_block_sizes);
+      this->column_block_indices.reinit(c_block_sizes);
+      this->sub_objects.reinit(nr, nc);
+      for (PetscInt i = 0; i < nr; ++i)
+        {
+          for (PetscInt j = 0; j < nc; ++j)
+            {
+              // TODO: MATNEST supports NULL blocks
+              this->sub_objects[i][j] = new BlockType(mats[i * nc + j]);
+            }
+        }
+
+      collect_sizes();
+    }
+
   } // namespace MPI
 } // namespace PETScWrappers
 

diff --git a/source/lac/petsc_parallel_sparse_matrix.cc b/source/lac/petsc_parallel_sparse_matrix.cc
index 7bbd637180cc6d5e6511998dc6d9f000599e5b61..6b6414b992c3b8f4d153328a2323c9d88cdec515 100644 (file)
--- a/source/lac/petsc_parallel_sparse_matrix.cc
+++ b/source/lac/petsc_parallel_sparse_matrix.cc
@@ -43,6 +43,9 @@ namespace PETScWrappers
       AssertThrow(ierr == 0, ExcPETScError(ierr));
     }
 
+    SparseMatrix::SparseMatrix(const Mat &A)
+      : MatrixBase(A)
+    {}
 
     SparseMatrix::~SparseMatrix()
     {

diff --git a/source/lac/petsc_sparse_matrix.cc b/source/lac/petsc_sparse_matrix.cc
index 194549925b485df272c68cdd89dfc4c2e026d06a..9d680bc0ac7082205c1b030abb3732ba97da3f83 100644 (file)
--- a/source/lac/petsc_sparse_matrix.cc
+++ b/source/lac/petsc_sparse_matrix.cc
@@ -35,7 +35,9 @@ namespace PETScWrappers
     AssertThrow(ierr == 0, ExcPETScError(ierr));
   }
 
-
+  SparseMatrix::SparseMatrix(const Mat &A)
+    : MatrixBase(A)
+  {}
 
   SparseMatrix::SparseMatrix(const size_type m,
                              const size_type n,

diff --git a/tests/petsc/block_matrices_01.cc b/tests/petsc/block_matrices_01.cc
index 193afc3b653e057f54d5bba132e34a2e3c18b0f8..d2a2fb9f81dfb4c2c1aef71ddbec324d40c74fbf 100644 (file)
--- a/tests/petsc/block_matrices_01.cc
+++ b/tests/petsc/block_matrices_01.cc
@@ -71,6 +71,26 @@ test()
   // Extract the PETSc MATNEST and use print from PETScWrappers::MatrixBase
   PETScWrappers::MatrixBase tmp(pbsm.petsc_matrix());
   tmp.print(deallog.get_file_stream());
+
+  // Extract the PETSc MATNEST and assign to a new BlockSparseMatrix
+  PETScWrappers::MPI::BlockSparseMatrix tmp2(pbsm.petsc_matrix());
+  Assert(tmp2.n_block_rows() == pbsm.n_block_rows(), ExcInternalError());
+  Assert(tmp2.n_block_cols() == pbsm.n_block_cols(), ExcInternalError());
+  Assert(tmp2.m() == pbsm.m(), ExcInternalError());
+  Assert(tmp2.n() == pbsm.n(), ExcInternalError());
+  for (unsigned int blr = 0; blr < 2; ++blr)
+    {
+      for (unsigned int blc = 0; blc < 2; ++blc)
+        {
+          Assert(tmp2.block(blr, blc).m() == pbsm.block(blr, blc).m(),
+                 ExcInternalError());
+          Assert(tmp2.block(blr, blc).n() == pbsm.block(blr, blc).n(),
+                 ExcInternalError());
+          Assert(tmp2.block(blr, blc).petsc_matrix() ==
+                   pbsm.block(blr, blc).petsc_matrix(),
+                 ExcInternalError());
+        }
+    }
 }