From e34ecec095746046bbf8a8ac3bd1517f098c96fb Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Thu, 10 Nov 2022 18:22:19 +0100 Subject: [PATCH] PETScWrappers:BlockSparseMatrix: create PETSc MATNEST --- .../deal.II/lac/petsc_block_sparse_matrix.h | 22 +++++++++++- .../lac/petsc_parallel_block_sparse_matrix.cc | 35 +++++++++++++++++++ tests/petsc/block_matrices_01.cc | 4 +++ tests/petsc/block_matrices_01.mpirun=1.output | 9 +++++ 4 files changed, 69 insertions(+), 1 deletion(-) diff --git a/include/deal.II/lac/petsc_block_sparse_matrix.h b/include/deal.II/lac/petsc_block_sparse_matrix.h index 913a694b87..8bb025586d 100644 --- a/include/deal.II/lac/petsc_block_sparse_matrix.h +++ b/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::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 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; }; diff --git a/source/lac/petsc_parallel_block_sparse_matrix.cc b/source/lac/petsc_parallel_block_sparse_matrix.cc index 6b82eae960..d744d72f2b 100644 --- a/source/lac/petsc_parallel_block_sparse_matrix.cc +++ b/source/lac/petsc_parallel_block_sparse_matrix.cc @@ -14,6 +14,7 @@ // --------------------------------------------------------------------- #include +#include #ifdef DEAL_II_WITH_PETSC @@ -31,6 +32,11 @@ namespace PETScWrappers return *this; } + BlockSparseMatrix::~BlockSparseMatrix() + { + PetscErrorCode ierr = destroy_matrix(petsc_nest_matrix); + AssertThrow(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 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 @@ -152,6 +176,17 @@ namespace PETScWrappers return block(0, 0).get_mpi_communicator(); } + BlockSparseMatrix::operator Mat() const + { + return petsc_nest_matrix; + } + + Mat & + BlockSparseMatrix::petsc_matrix() + { + return petsc_nest_matrix; + } + } // namespace MPI } // namespace PETScWrappers diff --git a/tests/petsc/block_matrices_01.cc b/tests/petsc/block_matrices_01.cc index 51dc0875dd..193afc3b65 100644 --- a/tests/petsc/block_matrices_01.cc +++ b/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()); } diff --git a/tests/petsc/block_matrices_01.mpirun=1.output b/tests/petsc/block_matrices_01.mpirun=1.output index a9c5e15408..ca11cbdd02 100644 --- a/tests/petsc/block_matrices_01.mpirun=1.output +++ b/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 -- 2.39.5