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 2fd3d37f0c1f26318a7bf6fe79f9048a583e0975..adbd203f27504a7ace68c54828dc098da781cfd1 100644 (file)
--- 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<size_type> &rows,
                const PetscScalar             new_diag_value = 0);
 
+    /**
+     * Same as clear_rows(), except that zeros also the columns.
+     */
+    void
+    clear_rows_columns(const std::vector<size_type> &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 b7dd5c688382dfb58c88b3b94b930158eb6372ac..83fecf13a565cd97466d110d56f0a6b3ae660312 100644 (file)
--- 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<size_type> &rows,
+                                 const PetscScalar             new_diag_value)
+  {
+    assert_is_compressed();
+
+    // now set all the entries of these rows
+    // to zero
+    const std::vector<PetscInt> 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 b86896cb64f2f5b4f6ca1bddee2d2cd369a0552c..2a4633cdc541f430973f6af95f333c0548a7d42d 100644 (file)
--- 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<types::global_dof_index> indices;
         std::vector<PetscScalar>             solution_values;
@@ -142,7 +143,10 @@ namespace MatrixTools
         // clear_rows() is a collective operation so we still have to call
         // it:
         std::vector<types::global_dof_index> 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 c7512144e0fbb39868b1ed85b73d5691d70821cf..431fab7531fecad46be609b7e56901425830ff2a 100644 (file)
--- 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 <deal.II/lac/petsc_sparse_matrix.h>
 #include <deal.II/lac/vector.h>
@@ -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<size_type>(&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 729838e63f09ec77846d62284c62057154ebc43f..8b29fc2b0beeffb6be4d424cc9870c228d44e44d 100644 (file)
--- 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 21b2e1348a597ad88f18c07e00842df77c940085..154139176fd0dd35d05eee1581aad38559cb1fdb 100644 (file)
--- 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 <deal.II/lac/petsc_sparse_matrix.h>
 #include <deal.II/lac/vector.h>
@@ -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<size_type>(&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 2000691fac310739f87466f06f3061ff2539c4df..35249b243ebe41ea67ded172535afea0ff070074 100644 (file)
--- 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