Allow creating matrices based on hash tables

clang-format

@@ -24,7 +24,7 @@ AllowShortLoopsOnASingleLine: false
 SortIncludes: false
 IndentCaseLabels: true
 ConstructorInitializerIndentWidth: 2
-AlwaysBreakAfterDefinitionReturnType: TopLevel
+AlwaysBreakAfterDefinitionReturnType: false
 AlwaysBreakTemplateDeclarations: true
 
 FixNamespaceComments: false

examples/systems_of_equations/systems_of_equations_ex7/systems_of_equations_ex7.C

@@ -563,6 +563,7 @@ int main (int argc, char ** argv)
 
   NonlinearImplicitSystem & system =
     equation_systems.add_system<NonlinearImplicitSystem> ("NonlinearElasticity");
+  system.prefer_hash_table_matrix_assembly(true);
 
   unsigned int u_var =
     system.add_variable("u",

include/numerics/petsc_matrix.h

@@ -76,14 +76,15 @@ class PetscMatrix final : public PetscMatrixBase<T>
 
   /**
    * Constructor.  Creates a PetscMatrix assuming you already have a
-   * valid Mat object.  In this case, m is NOT destroyed by the
+   * valid Mat object.  In this case, m may not be destroyed by the
    * PetscMatrix destructor when this object goes out of scope.  This
    * allows ownership of m to remain with the original creator, and to
    * simply provide additional functionality with the PetscMatrix.
    */
   explicit
   PetscMatrix (Mat m,
-               const Parallel::Communicator & comm_in);
+               const Parallel::Communicator & comm_in,
+               bool destroy_on_exit = false);
 
   /**
    * Constructor. Creates and initializes a PetscMatrix with the given
@@ -284,7 +285,48 @@ class PetscMatrix final : public PetscMatrixBase<T>
 
   virtual void scale(const T scale) override;
 
+#if PETSC_RELEASE_GREATER_EQUALS(3,23,0)
+  std::unique_ptr<PetscMatrix<T>> copy_from_hash();
+#endif
+
+  virtual bool supports_hash_table() const override;
+
+  virtual void reset_memory() override;
+
 protected:
+  /**
+   * Perform matrix initialization steps sans preallocation
+   * @param m The global number of rows
+   * @param n The global number of columns
+   * @param m_l The local number of rows
+   * @param n_l The local number of columns
+   * @param blocksize The matrix block size
+   */
+  void init_without_preallocation (numeric_index_type m,
+                                   numeric_index_type n,
+                                   numeric_index_type m_l,
+                                   numeric_index_type n_l,
+                                   numeric_index_type blocksize);
+
+  /*
+   * Performs matrix preallcation
+   * \param m_l The local number of rows.
+   * \param n_nz array containing the number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
+   * \param n_oz Array containing the number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
+   * \param blocksize Optional value indicating dense coupled blocks for systems with multiple variables all of the same    */
+  void preallocate(numeric_index_type m_l,
+                   const std::vector<numeric_index_type> & n_nz,
+                   const std::vector<numeric_index_type> & n_oz,
+                   numeric_index_type blocksize);
+
+  /**
+   * Finish up the initialization process. This method does a few things which include
+   * - Setting the option to make new nonzeroes an error (otherwise users will just have a silent
+       (often huge) performance penalty
+   * - Marking the matrix as initialized
+   * - Zeroing the matrix
+   */
+  void finish_initialization();
 
   /**
    * This function either creates or re-initializes a matrix called \p

include/numerics/petsc_matrix_base.h

@@ -84,14 +84,15 @@ class PetscMatrixBase : public SparseMatrix<T>
 
   /**
    * Constructor.  Creates a PetscMatrixBase assuming you already have a
-   * valid Mat object.  In this case, m is NOT destroyed by the
+   * valid Mat object.  In this case, m may not be destroyed by the
    * PetscMatrixBase destructor when this object goes out of scope.  This
    * allows ownership of m to remain with the original creator, and to
    * simply provide additional functionality with the PetscMatrixBase.
    */
   explicit
   PetscMatrixBase (Mat m,
-                   const Parallel::Communicator & comm_in);
+                   const Parallel::Communicator & comm_in,
+                   bool destroy_on_exit = false);
 
   /**
    * This class manages a C-style struct (Mat) manually, so we

include/numerics/sparse_matrix.h

@@ -583,6 +583,32 @@ class SparseMatrix : public ReferenceCountedObject<SparseMatrix<T>>,
    */
   virtual void scale(const T scale);
 
+  /**
+   * @returns Whether the matrix supports hash table assembly
+   */
+  virtual bool supports_hash_table() const { return false; }
+
+  /**
+   * @returns Whether a preference for hash table vs. preallocation has been set
+   */
+  bool assembly_preference_set() const { return _assembly_preference_set; }
+
+  /**
+   * Sets whether to use hash table assembly. This will error if the passed-in value is true and the
+   * matrix type does not support hash tables
+   */
+  void use_hash_table(bool use_hash);
+
+  /**
+   * @returns Whether this matrix is using hash table assembly
+   */
+  bool use_hash_table() const { return _use_hash_table; }
+
+  /**
+   * Reset the memory storage of the matrix
+   */
+  virtual void reset_memory() { libmesh_not_implemented(); }
+
 protected:
   /**
    * Protected implementation of the create_submatrix and reinit_submatrix
@@ -616,12 +642,33 @@ class SparseMatrix : public ReferenceCountedObject<SparseMatrix<T>>,
    * Flag indicating whether or not the matrix has been initialized.
    */
   bool _is_initialized;
+
+  /**
+   * Flag indicating whether the matrix is assembled using a hash table
+   */
+  bool _use_hash_table;
+
+  /**
+   * Whether a preference for hash table vs. preallocation has been set
+   */
+  bool _assembly_preference_set;
 };
 
 
 
 //-----------------------------------------------------------------------
 // SparseMatrix inline members
+template <typename T>
+void
+SparseMatrix<T>::use_hash_table(const bool use_hash)
+{
+  libmesh_error_msg_if(use_hash && !this->supports_hash_table(),
+                       "This matrix class does not support hash table assembly");
+  libmesh_error_msg_if(this->_assembly_preference_set,
+                       "The use_hash_table API should only be called once");
+  this->_use_hash_table = use_hash;
+  this->_assembly_preference_set = true;
+}
 
 // For SGI MIPSpro this implementation must occur after
 // the full specialization of the print() member.

include/systems/system.h

@@ -1922,6 +1922,10 @@ class System : public ReferenceCountedObject<System>,
    */
   SparseMatrix<Number> & get_matrix (std::string_view mat_name);
 
+  /**
+   * Sets whether to use hash table matrix assembly if the matrix sub-classes support it
+   */
+  void prefer_hash_table_matrix_assembly(bool preference);
 
 protected:
 
@@ -2293,6 +2297,16 @@ class System : public ReferenceCountedObject<System>,
    * Do we want to apply constraints while projecting vectors ?
    */
   bool project_with_constraints;
+
+  /**
+   * Whether to use hash table matrix assembly if the matrix sub-classes support it
+   */
+  bool _prefer_hash_table_matrix_assembly;
+
+  /**
+   * Whether any of our matrices require an initial sparsity pattern computation in order to determine preallocation
+   */
+  bool _require_sparsity_pattern;
 };
 
 
@@ -2685,6 +2699,14 @@ System::add_matrix (std::string_view mat_name,
   return mat;
 }
 
+inline void
+System::prefer_hash_table_matrix_assembly(const bool preference)
+{
+  libmesh_error_msg_if(
+      _matrices_initialized,
+      "System::prefer_hash_table_matrix_assembly() should be called before matrices are initialized");
+  _prefer_hash_table_matrix_assembly = preference;
+}
 
 } // namespace libMesh