From: Stefano Zampini Date: Sat, 26 Nov 2022 14:22:45 +0000 (+0300) Subject: PETScWrappers::MatrixBase implement clear_rows_columns X-Git-Tag: v9.5.0-rc1~635^2~2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=cb5892f7f8788f2966f3632579d66e96f0f1085f;p=dealii.git PETScWrappers::MatrixBase implement clear_rows_columns --- diff --git a/include/deal.II/lac/petsc_matrix_base.h b/include/deal.II/lac/petsc_matrix_base.h index 2fd3d37f0c..adbd203f27 100644 --- a/include/deal.II/lac/petsc_matrix_base.h +++ b/include/deal.II/lac/petsc_matrix_base.h @@ -571,6 +571,13 @@ namespace PETScWrappers clear_rows(const std::vector &rows, const PetscScalar new_diag_value = 0); + /** + * Same as clear_rows(), except that zeros also the columns. + */ + void + clear_rows_columns(const std::vector &rows, + const PetscScalar new_diag_value = 0); + /** * PETSc matrices store their own sparsity patterns. So, in analogy to our * own SparsityPattern class, this function compresses the sparsity diff --git a/source/lac/petsc_matrix_base.cc b/source/lac/petsc_matrix_base.cc index b7dd5c6883..83fecf13a5 100644 --- a/source/lac/petsc_matrix_base.cc +++ b/source/lac/petsc_matrix_base.cc @@ -168,6 +168,33 @@ namespace PETScWrappers ISDestroy(&index_set); } + void + MatrixBase::clear_rows_columns(const std::vector &rows, + const PetscScalar new_diag_value) + { + assert_is_compressed(); + + // now set all the entries of these rows + // to zero + const std::vector petsc_rows(rows.begin(), rows.end()); + + // call the functions. note that we have + // to call them even if #rows is empty, + // since this is a collective operation + IS index_set; + + ISCreateGeneral(get_mpi_communicator(), + rows.size(), + petsc_rows.data(), + PETSC_COPY_VALUES, + &index_set); + + const PetscErrorCode ierr = + MatZeroRowsColumnsIS(matrix, index_set, new_diag_value, nullptr, nullptr); + AssertThrow(ierr == 0, ExcPETScError(ierr)); + ISDestroy(&index_set); + } + PetscScalar diff --git a/source/numerics/matrix_tools_once.cc b/source/numerics/matrix_tools_once.cc index b86896cb64..2a4633cdc5 100644 --- a/source/numerics/matrix_tools_once.cc +++ b/source/numerics/matrix_tools_once.cc @@ -66,9 +66,6 @@ namespace MatrixTools PETScWrappers::VectorBase & right_hand_side, const bool eliminate_columns) { - (void)eliminate_columns; - Assert(eliminate_columns == false, ExcNotImplemented()); - Assert(matrix.n() == right_hand_side.size(), ExcDimensionMismatch(matrix.n(), right_hand_side.size())); Assert(matrix.n() == solution.size(), @@ -118,7 +115,11 @@ namespace MatrixTools // preserving it. this is different from // the case of deal.II sparse matrices // treated in the other functions. - matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry); + if (eliminate_columns) + matrix.clear_rows_columns(constrained_rows, + average_nonzero_diagonal_entry); + else + matrix.clear_rows(constrained_rows, average_nonzero_diagonal_entry); std::vector indices; std::vector solution_values; @@ -142,7 +143,10 @@ namespace MatrixTools // clear_rows() is a collective operation so we still have to call // it: std::vector constrained_rows; - matrix.clear_rows(constrained_rows, 1.); + if (eliminate_columns) + matrix.clear_rows_columns(constrained_rows, 1.); + else + matrix.clear_rows(constrained_rows, 1.); } // clean up @@ -159,6 +163,7 @@ namespace MatrixTools PETScWrappers::MPI::BlockVector & right_hand_side, const bool eliminate_columns) { + Assert(eliminate_columns == false, ExcNotImplemented()); Assert(matrix.n() == right_hand_side.size(), ExcDimensionMismatch(matrix.n(), right_hand_side.size())); Assert(matrix.n() == solution.size(), diff --git a/tests/petsc/67.cc b/tests/petsc/67.cc index c7512144e0..431fab7531 100644 --- a/tests/petsc/67.cc +++ b/tests/petsc/67.cc @@ -15,7 +15,8 @@ -// check PETScWrappers::MatrixBase::clear_rows () +// check PETScWrappers::MatrixBase::clear_rows () and +// PETScWrappers::MatrixBase::clear_rows_columns () #include #include @@ -97,6 +98,12 @@ test(PETScWrappers::MatrixBase &m) // may remove some, though) Assert(m.n_nonzero_elements() <= nnz, ExcInternalError()); + m.clear_rows_columns(std::vector(&rows[0], &rows[2])); + + deallog << m.frobenius_norm() << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + Assert(m.n_nonzero_elements() <= nnz, ExcInternalError()); deallog << "OK" << std::endl; } diff --git a/tests/petsc/67.output b/tests/petsc/67.output index 729838e63f..8b29fc2b0b 100644 --- a/tests/petsc/67.output +++ b/tests/petsc/67.output @@ -3,4 +3,6 @@ DEAL::7.15267e+09 7.15267e+09 DEAL::32 32 DEAL::6.71272e+09 6.71272e+09 DEAL::27 32 +DEAL::6.04477e+09 +DEAL::27 32 DEAL::OK diff --git a/tests/petsc/69.cc b/tests/petsc/69.cc index 21b2e1348a..154139176f 100644 --- a/tests/petsc/69.cc +++ b/tests/petsc/69.cc @@ -15,7 +15,8 @@ -// check PETScWrappers::MatrixBase::clear_rows () with used second argument +// check PETScWrappers::MatrixBase::clear_rows () and +// PETScWrappers::MatrixBase::clear_rows_columns () with used second argument #include #include @@ -101,6 +102,12 @@ test(PETScWrappers::MatrixBase &m) // may remove some, though) AssertThrow(m.n_nonzero_elements() <= nnz, ExcInternalError()); + m.clear_rows_columns(std::vector(&rows[0], &rows[2]), rnd); + + deallog << m.frobenius_norm() << std::endl; + deallog << m.n_nonzero_elements() << ' ' << nnz << std::endl; + + AssertThrow(m.n_nonzero_elements() <= nnz, ExcInternalError()); deallog << "OK" << std::endl; } diff --git a/tests/petsc/69.output b/tests/petsc/69.output index 2000691fac..35249b243e 100644 --- a/tests/petsc/69.output +++ b/tests/petsc/69.output @@ -3,4 +3,6 @@ DEAL::7.15267e+09 7.15267e+09 DEAL::32 32 DEAL::6.96580e+09 6.96580e+09 DEAL::29 32 +DEAL::6.32463e+09 +DEAL::29 32 DEAL::OK