PETScWrappers::MatrixBase implement clear_rows_columns

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

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

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(),

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;
 }
 

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

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;
 }
 

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