Merge pull request #14592 from luca-heltai/petsc-block-matrix-to-nestmat

PETScWrappers:BlockSparseMatrix: create PETSc MATNEST

doc/news/changes/minor/20221220StefanoZampini-b

@@ -0,0 +1,3 @@
+New: PETScWrappers:BlockSparseMatrix now is also a PETSc MATNEST type.
+<br>
+(Stefano Zampini, 2022/12/20)

include/deal.II/lac/petsc_block_sparse_matrix.h

@@ -105,7 +105,7 @@ namespace PETScWrappers
       /**
        * Destructor.
        */
-      ~BlockSparseMatrix() override = default;
+      ~BlockSparseMatrix() override;
 
       /**
        * Pseudo copy operator only copying empty objects. The sizes of the
@@ -271,6 +271,26 @@ namespace PETScWrappers
        * protected.
        */
       using BlockMatrixBase<SparseMatrix>::clear;
+
+      /**
+       * Conversion operator to gain access to the underlying PETSc type. If you
+       * do this, you cut this class off some information it may need, so this
+       * conversion operator should only be used if you know what you do. In
+       * particular, it should only be used for read-only operations into the
+       * matrix.
+       */
+      operator const 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();
+
+    private:
+      Mat petsc_nest_matrix = nullptr;
     };
 
 

source/lac/petsc_parallel_block_sparse_matrix.cc

@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/lac/petsc_block_sparse_matrix.h>
+#include <deal.II/lac/petsc_compatibility.h>
 
 #ifdef DEAL_II_WITH_PETSC
 
@@ -31,6 +32,11 @@ namespace PETScWrappers
       return *this;
     }
 
+    BlockSparseMatrix::~BlockSparseMatrix()
+    {
+      PetscErrorCode ierr = destroy_matrix(petsc_nest_matrix);
+      AssertNothrow(ierr == 0, ExcPETScError(ierr));
+    }
 
 #  ifndef DOXYGEN
     void
@@ -111,6 +117,24 @@ namespace PETScWrappers
     BlockSparseMatrix::collect_sizes()
     {
       BaseClass::collect_sizes();
+
+      auto m = this->n_block_cols();
+      auto n = this->n_block_cols();
+
+      PetscErrorCode ierr = destroy_matrix(petsc_nest_matrix);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
+      std::vector<Mat> psub_objects(m * n);
+      for (unsigned int r = 0; r < m; r++)
+        for (unsigned int c = 0; c < n; c++)
+          psub_objects[r * n + c] = this->sub_objects[r][c]->petsc_matrix();
+      ierr = MatCreateNest(get_mpi_communicator(),
+                           m,
+                           NULL,
+                           n,
+                           NULL,
+                           psub_objects.data(),
+                           &petsc_nest_matrix);
+      AssertThrow(ierr == 0, ExcPETScError(ierr));
     }
 
     std::vector<IndexSet>
@@ -152,6 +176,17 @@ namespace PETScWrappers
       return block(0, 0).get_mpi_communicator();
     }
 
+    BlockSparseMatrix::operator const Mat &() const
+    {
+      return petsc_nest_matrix;
+    }
+
+    Mat &
+    BlockSparseMatrix::petsc_matrix()
+    {
+      return petsc_nest_matrix;
+    }
+
   } // namespace MPI
 } // namespace PETScWrappers
 

tests/petsc/block_matrices_01.cc

@@ -67,6 +67,10 @@ test()
   pbsm.reinit(rows, cols, bdsp, MPI_COMM_WORLD);
   deallog << "nonzeros BlockSparseMatrix: " << pbsm.n_nonzero_elements()
           << std::endl;
+
+  // Extract the PETSc MATNEST and use print from PETScWrappers::MatrixBase
+  PETScWrappers::MatrixBase tmp(pbsm.petsc_matrix());
+  tmp.print(deallog.get_file_stream());
 }
 
 

tests/petsc/block_matrices_01.mpirun=1.output

@@ -1,3 +1,12 @@
 
 DEAL::nonzeros BlockSparsityPattern: 9
 DEAL::nonzeros BlockSparseMatrix: 9
+(0,0) 0.00000
+(1,0) 0.00000
+(1,1) 0.00000
+(2,2) 0.00000
+(3,2) 0.00000
+(3,3) 0.00000
+(4,2) 0.00000
+(4,3) 0.00000
+(4,4) 0.00000