From 3149b0647edd6feb013e987c4e9b59ef6f01f965 Mon Sep 17 00:00:00 2001 From: Stefano Zampini Date: Mon, 14 Nov 2022 21:33:39 +0100 Subject: [PATCH] PETScWrappers:BlockSparseMatrix: constructor from PETSc Mat --- .../deal.II/lac/petsc_block_sparse_matrix.h | 17 +++++++ include/deal.II/lac/petsc_matrix_base.h | 7 +++ include/deal.II/lac/petsc_sparse_matrix.h | 12 +++++ source/lac/petsc_matrix_base.cc | 13 ++++- .../lac/petsc_parallel_block_sparse_matrix.cc | 50 +++++++++++++++++++ source/lac/petsc_parallel_sparse_matrix.cc | 3 ++ source/lac/petsc_sparse_matrix.cc | 4 +- tests/petsc/block_matrices_01.cc | 20 ++++++++ 8 files changed, 123 insertions(+), 3 deletions(-) diff --git a/include/deal.II/lac/petsc_block_sparse_matrix.h b/include/deal.II/lac/petsc_block_sparse_matrix.h index f62d718b2a..626974e70a 100644 --- 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 2fd3d37f0c..366d4172ee 100644 --- 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 08d7734571..a3ced65c7c 100644 --- 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 b7dd5c6883..c7363ba9c7 100644 --- 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(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 a290cfb602..ae50044645 100644 --- 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(A), + MATNEST, + &isnest); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + std::vector 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 r_block_sizes(nr, 0); + std::vector 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 7bbd637180..6b6414b992 100644 --- 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 194549925b..9d680bc0ac 100644 --- 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 193afc3b65..d2a2fb9f81 100644 --- 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()); + } + } } -- 2.39.5