Merge pull request #14684 from luca-heltai/stefanozampini/petsc-zeror…

…owscolumns

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 the function also zeros the columns.
+     */
+    void
+    clear_rows_columns(const std::vector<size_type> &row_and_column_indices,
+                       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

include/deal.II/numerics/matrix_tools.h

@@ -345,35 +345,7 @@ namespace MatrixTools
    * described in the general documentation of this namespace. This function
    * works on the classes that are used to wrap PETSc objects.
    *
-   * <b>Important:</b> This function is not very efficient: it needs to
-   * alternatingly read and write into the matrix, a situation that PETSc does
-   * not handle well. In addition, we only get rid of rows corresponding to
-   * boundary nodes, but the corresponding case of deleting the respective
-   * columns (i.e. if @p eliminate_columns is @p true) is not presently
-   * implemented, and probably will never because it is too expensive without
-   * direct access to the PETSc data structures. (This leads to the situation
-   * where the action indicated by the default value of the last argument is
-   * actually not implemented; that argument has <code>true</code> as its
-   * default value to stay consistent with the other functions of same name in
-   * this namespace.)
-   *
    * This function is used in step-17 and step-18.
-   *
-   * @note If the matrix is stored in parallel across multiple processors
-   * using MPI, this function only touches rows that are locally stored and
-   * simply ignores all other rows. In other words, each processor is
-   * responsible for its own rows, and the @p boundary_values argument needs
-   * to contain all locally owned rows of the matrix that you want to have
-   * treated. (But it can also contain entries for degrees of freedom not
-   * owned locally; these will simply be ignored.) Further, in the context of
-   * parallel computations, you will get into trouble if you treat a row while
-   * other processors still have pending writes or additions into the same
-   * row. In other words, if another processor still wants to add something to
-   * an element of a row and you call this function to zero out the row, then
-   * the next time you call compress() may add the remote value to the zero
-   * you just created. Consequently, you will want to call compress() after
-   * you made the last modifications to a matrix and before starting to clear
-   * rows.
    */
   void
   apply_boundary_values(

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