Merge pull request #2169 from lindsayad/expand-petsc-interface

Expand PetscMatrix interface

Author: lindsayad

include/numerics/petsc_matrix.h

@@ -151,6 +151,11 @@ class PetscMatrix final : public SparseMatrix<T>
    */
   void update_preallocation_and_zero();
 
+  /**
+   * Reset matrix to use the original nonzero pattern provided by users
+   */
+  void reset_preallocation();
+
   virtual void clear () override;
 
   virtual void zero () override;
@@ -163,8 +168,25 @@ class PetscMatrix final : public SparseMatrix<T>
 
   virtual numeric_index_type m () const override;
 
+  /**
+   * Get the number of rows owned by this process
+   */
+  numeric_index_type local_m () const;
+
   virtual numeric_index_type n () const override;
 
+  /**
+   * Get the number of columns owned by this process
+   */
+  numeric_index_type local_n () const;
+
+  /**
+   * Get the number of rows and columns owned by this process
+   * @param i Row size
+   * @param j Column size
+   */
+  void get_local_size (numeric_index_type & m, numeric_index_type & n) const;
+
   virtual numeric_index_type row_start () const override;
 
   virtual numeric_index_type row_stop () const override;

src/numerics/petsc_matrix.C

@@ -172,11 +172,11 @@ void PetscMatrix<T>::init (const numeric_index_type m_in,
 
       ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
       LIBMESH_CHKERR(ierr);
-      ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, PETSC_NULL);
+      ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, NULL);
       LIBMESH_CHKERR(ierr);
       ierr = MatMPIBAIJSetPreallocation(_mat, blocksize,
-                                        n_nz/blocksize, PETSC_NULL,
-                                        n_oz/blocksize, PETSC_NULL);
+                                        n_nz/blocksize, NULL,
+                                        n_oz/blocksize, NULL);
       LIBMESH_CHKERR(ierr);
     }
   else
@@ -186,17 +186,17 @@ void PetscMatrix<T>::init (const numeric_index_type m_in,
         case AIJ:
           ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
           LIBMESH_CHKERR(ierr);
-          ierr = MatSeqAIJSetPreallocation(_mat, n_nz, PETSC_NULL);
+          ierr = MatSeqAIJSetPreallocation(_mat, n_nz, NULL);
           LIBMESH_CHKERR(ierr);
-          ierr = MatMPIAIJSetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
+          ierr = MatMPIAIJSetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
           LIBMESH_CHKERR(ierr);
           break;
 
         case HYPRE:
 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
           ierr = MatSetType(_mat, MATHYPRE);
           LIBMESH_CHKERR(ierr);
-          ierr = MatHYPRESetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
+          ierr = MatHYPRESetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
           LIBMESH_CHKERR(ierr);
 #else
           libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
@@ -470,6 +470,20 @@ void PetscMatrix<T>::update_preallocation_and_zero ()
   libmesh_not_implemented();
 }
 
+template <typename T>
+void PetscMatrix<T>::reset_preallocation()
+{
+#if !PETSC_VERSION_LESS_THAN(3,8,0)
+  libmesh_assert (this->initialized());
+
+  auto ierr = MatResetPreallocation(_mat);
+  LIBMESH_CHKERR(ierr);
+#else
+  libmesh_warning("Your version of PETSc doesn't support resetting of "
+                  "preallocation, so we will use your most recent sparsity "
+                  "pattern. This may result in a degradation of performance\n");
+#endif
+}
 
 template <typename T>
 void PetscMatrix<T>::zero ()
@@ -506,7 +520,7 @@ void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_v
     ierr = MatZeroRows(_mat, rows.size(),
                        numeric_petsc_cast(rows.data()), diag_value);
   else
-    ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value);
+    ierr = MatZeroRows(_mat, 0, NULL, diag_value);
 #else
   // As of petsc-dev at the time of 3.1.0, MatZeroRows now takes two additional
   // optional arguments.  The optional arguments (x,b) can be used to specify the
@@ -515,10 +529,10 @@ void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_v
   if (!rows.empty())
     ierr = MatZeroRows(_mat, cast_int<PetscInt>(rows.size()),
                        numeric_petsc_cast(rows.data()), diag_value,
-                       PETSC_NULL, PETSC_NULL);
+                       NULL, NULL);
   else
-    ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value, PETSC_NULL,
-                       PETSC_NULL);
+    ierr = MatZeroRows(_mat, 0, NULL, diag_value, NULL,
+                       NULL);
 #endif
 
   LIBMESH_CHKERR(ierr);
@@ -931,7 +945,7 @@ void PetscMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
   PetscErrorCode ierr;
 #if PETSC_VERSION_LESS_THAN(3,0,0)
   if (&petsc_dest == this)
-    ierr = MatTranspose(_mat,PETSC_NULL);
+    ierr = MatTranspose(_mat,NULL);
   else
     ierr = MatTranspose(_mat,&petsc_dest._mat);
   LIBMESH_CHKERR(ierr);
@@ -1002,7 +1016,18 @@ numeric_index_type PetscMatrix<T>::m () const
   return static_cast<numeric_index_type>(petsc_m);
 }
 
+template <typename T>
+numeric_index_type PetscMatrix<T>::local_m () const
+{
+  libmesh_assert (this->initialized());
+
+  PetscInt m = 0;
 
+  auto ierr = MatGetLocalSize (_mat, &m, NULL);
+  LIBMESH_CHKERR(ierr);
+
+  return static_cast<numeric_index_type>(m);
+}
 
 template <typename T>
 numeric_index_type PetscMatrix<T>::n () const
@@ -1018,7 +1043,33 @@ numeric_index_type PetscMatrix<T>::n () const
   return static_cast<numeric_index_type>(petsc_n);
 }
 
+template <typename T>
+numeric_index_type PetscMatrix<T>::local_n () const
+{
+  libmesh_assert (this->initialized());
+
+  PetscInt n = 0;
+
+  auto ierr = MatGetLocalSize (_mat, NULL, &n);
+  LIBMESH_CHKERR(ierr);
+
+  return static_cast<numeric_index_type>(n);
+}
 
+template <typename T>
+void PetscMatrix<T>::get_local_size (numeric_index_type & m,
+                                     numeric_index_type & n) const
+{
+  libmesh_assert (this->initialized());
+
+  PetscInt petsc_m = 0, petsc_n = 0;
+
+  auto ierr = MatGetLocalSize (_mat, &petsc_m, &petsc_n);
+  LIBMESH_CHKERR(ierr);
+
+  m = static_cast<numeric_index_type>(petsc_m);
+  n = static_cast<numeric_index_type>(petsc_n);
+}
 
 template <typename T>
 numeric_index_type PetscMatrix<T>::row_start () const