Miscennaleous Example 17 - Demonstrating mix of preallocated/hash-table matrix
+// assemblies
\author Alexander D. Lindsay \date 2024
+
+// Basic include files needed for the mesh functionality.
+#include "libmesh/libmesh.h"
+#include "libmesh/mesh.h"
+#include "libmesh/mesh_generation.h"
+#include "libmesh/linear_implicit_system.h"
+#include "libmesh/equation_systems.h"
+
+// Define the Finite Element object.
+#include "libmesh/fe.h"
+
+// Define Gauss quadrature rules.
+#include "libmesh/quadrature_gauss.h"
+
+// Define useful datatypes for finite element
+// matrix and vector components.
+#include "libmesh/numeric_vector.h"
+#include "libmesh/sparse_matrix.h"
+#include "libmesh/dense_matrix.h"
+#include "libmesh/dense_vector.h"
+#include "libmesh/elem.h"
+#include "libmesh/enum_solver_package.h"
+
+// Define the DofMap, which handles degree of freedom
+// indexing.
+#include "libmesh/dof_map.h"
+
+#ifdef LIBMESH_HAVE_PETSC
+// include PETSc headers
+#include "libmesh/petsc_matrix.h"
+#include "libmesh/petsc_vector.h"
+#include "petscksp.h"
+#endif
+
+#include
+
+// Bring in everything from the libMesh namespace
+using namespace libMesh;
+
+void assemble_poisson(EquationSystems & es,
+ const std::string & system_name);
+
+// Function prototype for the exact solution.
+Real exact_solution (const Real x,
+ const Real y,
+ const Real z = 0.);
+
+int main (int argc, char ** argv)
+{
+ // Initialize libraries, like in example 2.
+ LibMeshInit init (argc, argv);
+
+ // Brief message to the user regarding the program name
+ // and command line arguments.
+ libMesh::out << "Running " << argv[0];
+
+ for (int i=1; i("Poisson");
+
+ // Adds the variable "u" to "Poisson". "u"
+ // will be approximated using second-order approximation.
+ system.add_variable("u", FIRST);
+
+ // Give the system a pointer to the matrix assembly
+ // function. This will be called when needed by the
+ // library.
+ system.attach_assemble_function(assemble_poisson);
+
+ // Add the preconditioner matrix
+ system.add_matrix("preconditioner");
+#ifdef LIBMESH_HAVE_PETSC
+#if PETSC_RELEASE_GREATER_EQUALS(3, 19, 0)
+ system.get_matrix("preconditioner").use_hash_table(true);
+#endif
+#endif
+
+ // Initialize the data structures for the equation system.
+ equation_systems.init();
+
+ // Prints information about the system to the screen.
+ equation_systems.print_info();
+
+ // assemble the operators and RHS
+ system.assemble();
+
+#ifdef LIBMESH_HAVE_PETSC
+ auto & sys_matrix = cast_ref &>(system.get_system_matrix());
+ auto & pre_matrix = cast_ref &>(system.get_matrix("preconditioner"));
+ LibmeshPetscCall2(system.comm(), PetscOptionsSetValue(NULL, "-ksp_monitor", NULL));
+
+ auto solve = [&sys_matrix, &pre_matrix, &system]() {
+ // Make sure our matrices are closed
+ sys_matrix.close();
+ pre_matrix.close();
+ system.rhs->close();
+ KSP ksp;
+ LibmeshPetscCall2(system.comm(), KSPCreate(system.comm().get(), &ksp));
+ LibmeshPetscCall2(system.comm(), KSPSetOperators(ksp, sys_matrix.mat(), pre_matrix.mat()));
+ LibmeshPetscCall2(system.comm(), KSPSetFromOptions(ksp));
+ LibmeshPetscCall2(system.comm(),
+ KSPSolve(ksp,
+ cast_ptr *>(system.rhs)->vec(),
+ cast_ptr *>(system.solution.get())->vec()));
+ };
+
+ // solve
+ solve();
+
+ // MatResetHash added in PETSc version 3.23
+#if !PETSC_VERSION_LESS_THAN(3, 23, 0)
+ // reset the memory
+ // sys_matrix.reset_memory(); # See https://gitlab.com/petsc/petsc/-/merge_requests/8063
+ pre_matrix.reset_memory();
+ // zero
+ sys_matrix.zero();
+ pre_matrix.zero();
+ system.rhs->zero();
+ system.solution->zero();
+ system.update();
+ // re-assemble
+ system.assemble();
+ // resolve
+ solve();
+#endif
+#endif
+
+ // All done.
+ return 0;
+}
+
+
+
+// We now define the matrix assembly function for the
+// Poisson system. We need to first compute element
+// matrices and right-hand sides, and then take into
+// account the boundary conditions, which will be handled
+// via a penalty method.
+void assemble_poisson(EquationSystems & es,
+ const std::string & libmesh_dbg_var(system_name))
+{
+
+ // It is a good idea to make sure we are assembling
+ // the proper system.
+ libmesh_assert_equal_to (system_name, "Poisson");
+
+ // Get a constant reference to the mesh object.
+ const MeshBase & mesh = es.get_mesh();
+
+ // The dimension that we are running
+ const unsigned int dim = mesh.mesh_dimension();
+
+ // Get a reference to the LinearImplicitSystem we are solving
+ LinearImplicitSystem & system = es.get_system ("Poisson");
+
+ // A reference to the DofMap object for this system. The DofMap
+ // object handles the index translation from node and element numbers
+ // to degree of freedom numbers. We will talk more about the DofMap
+ // in future examples.
+ const DofMap & dof_map = system.get_dof_map();
+
+ // Get a constant reference to the Finite Element type
+ // for the first (and only) variable in the system.
+ FEType fe_type = dof_map.variable_type(0);
+
+ // Build a Finite Element object of the specified type. Since the
+ // FEBase::build() member dynamically creates memory we will
+ // store the object as a std::unique_ptr. This can be thought
+ // of as a pointer that will clean up after itself. Introduction Example 4
+ // describes some advantages of std::unique_ptr's in the context of
+ // quadrature rules.
+ std::unique_ptr fe (FEBase::build(dim, fe_type));
+
+ // A 5th order Gauss quadrature rule for numerical integration.
+ QGauss qrule (dim, FIFTH);
+
+ // Tell the finite element object to use our quadrature rule.
+ fe->attach_quadrature_rule (&qrule);
+
+ // Declare a special finite element object for
+ // boundary integration.
+ std::unique_ptr fe_face (FEBase::build(dim, fe_type));
+
+ // Boundary integration requires one quadrature rule,
+ // with dimensionality one less than the dimensionality
+ // of the element.
+ QGauss qface(dim-1, FIFTH);
+
+ // Tell the finite element object to use our
+ // quadrature rule.
+ fe_face->attach_quadrature_rule (&qface);
+
+ // Here we define some references to cell-specific data that
+ // will be used to assemble the linear system.
+ //
+ // The element Jacobian * quadrature weight at each integration point.
+ const std::vector & JxW = fe->get_JxW();
+
+ // The physical XY locations of the quadrature points on the element.
+ // These might be useful for evaluating spatially varying material
+ // properties at the quadrature points.
+ const std::vector & q_point = fe->get_xyz();
+
+ // The element shape functions evaluated at the quadrature points.
+ const std::vector> & phi = fe->get_phi();
+
+ // The element shape function gradients evaluated at the quadrature
+ // points.
+ const std::vector> & dphi = fe->get_dphi();
+
+ // Define data structures to contain the element matrix
+ // and right-hand-side vector contribution. Following
+ // basic finite element terminology we will denote these
+ // "Ke" and "Fe". These datatypes are templated on
+ // Number, which allows the same code to work for real
+ // or complex numbers.
+ DenseMatrix Ke;
+ DenseVector Fe;
+
+ // This vector will hold the degree of freedom indices for
+ // the element. These define where in the global system
+ // the element degrees of freedom get mapped.
+ std::vector dof_indices;
+
+ // The global system matrix
+ SparseMatrix & matrix = system.get_system_matrix();
+ // The preconditioning matrix
+ auto & pre_matrix = system.get_matrix("preconditioner");
+
+ // Now we will loop over all the elements in the mesh.
+ // We will compute the element matrix and right-hand-side
+ // contribution.
+ //
+ // Element ranges are a nice way to iterate through all the
+ // elements, or all the elements that have some property. The
+ // range will iterate from the first to the last element on
+ // the local processor.
+ // It is smart to make this one const so that we don't accidentally
+ // mess it up! In case users later modify this program to include
+ // refinement, we will be safe and will only consider the active
+ // elements; hence we use a variant of the
+ // active_local_element_ptr_range.
+ for (const auto & elem : mesh.active_local_element_ptr_range())
+ {
+ // Get the degree of freedom indices for the
+ // current element. These define where in the global
+ // matrix and right-hand-side this element will
+ // contribute to.
+ dof_map.dof_indices (elem, dof_indices);
+
+ // Cache the number of degrees of freedom on this element, for
+ // use as a loop bound later. We use cast_int to explicitly
+ // convert from size() (which may be 64-bit) to unsigned int
+ // (which may be 32-bit but which is definitely enough to count
+ // *local* degrees of freedom.
+ const unsigned int n_dofs =
+ cast_int(dof_indices.size());
+
+ // Compute the element-specific data for the current
+ // element. This involves computing the location of the
+ // quadrature points (q_point) and the shape functions
+ // (phi, dphi) for the current element.
+ fe->reinit (elem);
+
+ // With one variable, we should have the same number of degrees
+ // of freedom as shape functions.
+ libmesh_assert_equal_to (n_dofs, phi.size());
+
+ // Zero the element matrix and right-hand side before
+ // summing them. We use the resize member here because
+ // the number of degrees of freedom might have changed from
+ // the last element. Note that this will be the case if the
+ // element type is different (i.e. the last element was a
+ // triangle, now we are on a quadrilateral).
+
+ // The DenseMatrix::resize() and the DenseVector::resize()
+ // members will automatically zero out the matrix and vector.
+ Ke.resize (n_dofs, n_dofs);
+
+ Fe.resize (n_dofs);
+
+ // Now loop over the quadrature points. This handles
+ // the numeric integration.
+ for (unsigned int qp=0; qpside_index_range())
+ if (elem->neighbor_ptr(side) == nullptr)
+ {
+ // The value of the shape functions at the quadrature
+ // points.
+ const std::vector> & phi_face = fe_face->get_phi();
+
+ // The Jacobian * Quadrature Weight at the quadrature
+ // points on the face.
+ const std::vector & JxW_face = fe_face->get_JxW();
+
+ // The XYZ locations (in physical space) of the
+ // quadrature points on the face. This is where
+ // we will interpolate the boundary value function.
+ const std::vector & qface_point = fe_face->get_xyz();
+
+ // Compute the shape function values on the element
+ // face.
+ fe_face->reinit(elem, side);
+
+ // Some shape functions will be 0 on the face, but for
+ // ease of indexing and generality of code we loop over
+ // them anyway
+ libmesh_assert_equal_to (n_dofs, phi_face.size());
+
+ // Loop over the face quadrature points for integration.
+ for (unsigned int qp=0; qpadd_vector (Fe, dof_indices);
+ }
+
+ // All done!
+}
diff --git a/examples/miscellaneous/miscellaneous_ex17/run.sh b/examples/miscellaneous/miscellaneous_ex17/run.sh
new file mode 100755
index 00000000000..1af988300cc
--- /dev/null
+++ b/examples/miscellaneous/miscellaneous_ex17/run.sh
@@ -0,0 +1,9 @@
+#!/bin/sh
+
+#set -x
+
+. "$LIBMESH_DIR"/examples/run_common.sh
+
+example_name=miscellaneous_ex17
+
+run_example "$example_name"
diff --git a/examples/systems_of_equations/systems_of_equations_ex7/systems_of_equations_ex7.C b/examples/systems_of_equations/systems_of_equations_ex7/systems_of_equations_ex7.C
index 436b9d92433..19fc37a7fde 100644
--- a/examples/systems_of_equations/systems_of_equations_ex7/systems_of_equations_ex7.C
+++ b/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 ("NonlinearElasticity");
+ system.prefer_hash_table_matrix_assembly(true);
unsigned int u_var =
system.add_variable("u",
diff --git a/include/numerics/petsc_matrix.h b/include/numerics/petsc_matrix.h
index e7688eb5853..7d6f7a4fbee 100644
--- a/include/numerics/petsc_matrix.h
+++ b/include/numerics/petsc_matrix.h
@@ -36,7 +36,7 @@
#include
#endif
-
+class PetscMatrixTest;
namespace libMesh
{
@@ -76,14 +76,15 @@ class PetscMatrix final : public PetscMatrixBase
/**
* 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,52 @@ class PetscMatrix final : public PetscMatrixBase
virtual void scale(const T scale) override;
+#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
+ /**
+ * Creates a copy of the current hash table matrix and then performs assembly. This is very useful
+ * in cases where you are not done filling this matrix but want to be able to read the current
+ * state of it
+ */
+ std::unique_ptr> 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 & n_nz,
+ const std::vector & 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
+ */
+ void finish_initialization();
/**
* This function either creates or re-initializes a matrix called \p
@@ -313,6 +359,8 @@ class PetscMatrix final : public PetscMatrixBase
#else
mutable Threads::spin_mutex _petsc_matrix_mutex;
#endif
+
+ friend class ::PetscMatrixTest;
};
} // namespace libMesh
diff --git a/include/numerics/petsc_matrix_base.h b/include/numerics/petsc_matrix_base.h
index cda01c875eb..61124433132 100644
--- a/include/numerics/petsc_matrix_base.h
+++ b/include/numerics/petsc_matrix_base.h
@@ -84,14 +84,15 @@ class PetscMatrixBase : public SparseMatrix
/**
* 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
diff --git a/include/numerics/sparse_matrix.h b/include/numerics/sparse_matrix.h
index 840d1d5e829..1f9a58b07d0 100644
--- a/include/numerics/sparse_matrix.h
+++ b/include/numerics/sparse_matrix.h
@@ -583,6 +583,33 @@ class SparseMatrix : public ReferenceCountedObject>,
*/
virtual void scale(const T scale);
+ /**
+ * @returns Whether the matrix supports hash table assembly
+ */
+ virtual bool supports_hash_table() const { return false; }
+
+ /**
+ * 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. Hash table or hash map assembly means storing maps
+ * from i-j locations in the matrix to values. Because it is a hash map as opposed to a contiguous
+ * array of data, no preallocation is required to use it
+ */
+ void use_hash_table(bool use_hash);
+
+ /**
+ * @returns Whether this matrix is using hash table assembly. Hash table or hash map assembly
+ * means storing maps from i-j locations in the matrix to values. Because it is a hash map as
+ * opposed to a contiguous array of data, no preallocation is required to use it
+ */
+ bool use_hash_table() const { return _use_hash_table; }
+
+ /**
+ * Reset the memory storage of the matrix. Unlike \p clear(), this does not destroy the matrix but
+ * rather will reset the matrix to use the original preallocation or when using hash table matrix
+ * assembly (see \p use_hash_table()) will reset (clear) the hash table used for assembly
+ */
+ virtual void reset_memory() { libmesh_not_implemented(); }
+
protected:
/**
* Protected implementation of the create_submatrix and reinit_submatrix
@@ -616,12 +643,25 @@ class SparseMatrix : public ReferenceCountedObject>,
* 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;
};
//-----------------------------------------------------------------------
// SparseMatrix inline members
+template
+void
+SparseMatrix::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");
+ this->_use_hash_table = use_hash;
+}
// For SGI MIPSpro this implementation must occur after
// the full specialization of the print() member.
diff --git a/include/systems/system.h b/include/systems/system.h
index e7fa2442119..26215fdd70e 100644
--- a/include/systems/system.h
+++ b/include/systems/system.h
@@ -1922,6 +1922,10 @@ class System : public ReferenceCountedObject,
*/
SparseMatrix & 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,
* 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
diff --git a/src/numerics/petsc_matrix.C b/src/numerics/petsc_matrix.C
index f38801a5dfd..a5dd9ff6a37 100644
--- a/src/numerics/petsc_matrix.C
+++ b/src/numerics/petsc_matrix.C
@@ -94,8 +94,9 @@ PetscMatrix::PetscMatrix(const Parallel::Communicator & comm_in) :
// for destroying it
template
PetscMatrix::PetscMatrix(Mat mat_in,
- const Parallel::Communicator & comm_in) :
- PetscMatrixBase(mat_in, comm_in)
+ const Parallel::Communicator & comm_in,
+ const bool destroy_on_exit) :
+ PetscMatrixBase(mat_in, comm_in, destroy_on_exit)
{
MatType mat_type;
LibmeshPetscCall(MatGetType(mat_in, &mat_type));
@@ -129,16 +130,13 @@ PetscMatrix::PetscMatrix(const Parallel::Communicator & comm_in,
template
PetscMatrix::~PetscMatrix() = default;
-
-
template
-void PetscMatrix::init (const numeric_index_type m_in,
- const numeric_index_type n_in,
- const numeric_index_type m_l,
- const numeric_index_type n_l,
- const numeric_index_type nnz,
- const numeric_index_type noz,
- const numeric_index_type blocksize_in)
+void
+PetscMatrix::init_without_preallocation (const numeric_index_type m_in,
+ const numeric_index_type n_in,
+ const numeric_index_type m_l,
+ const numeric_index_type n_l,
+ const numeric_index_type blocksize_in)
{
// So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
libmesh_ignore(blocksize_in);
@@ -147,20 +145,16 @@ void PetscMatrix::init (const numeric_index_type m_in,
if (this->initialized())
this->clear();
- this->_is_initialized = true;
-
-
PetscInt m_global = static_cast(m_in);
PetscInt n_global = static_cast(n_in);
PetscInt m_local = static_cast(m_l);
PetscInt n_local = static_cast(n_l);
- PetscInt n_nz = static_cast(nnz);
- PetscInt n_oz = static_cast(noz);
LibmeshPetscCall(MatCreate(this->comm().get(), &this->_mat));
LibmeshPetscCall(MatSetSizes(this->_mat, m_local, n_local, m_global, n_global));
PetscInt blocksize = static_cast(blocksize_in);
LibmeshPetscCall(MatSetBlockSize(this->_mat,blocksize));
+ LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
if (blocksize > 1)
@@ -179,12 +173,7 @@ void PetscMatrix::init (const numeric_index_type m_in,
// MatSetFromOptions needs to happen before Preallocation routines
// since MatSetFromOptions can change matrix type and remove incompatible
// preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
LibmeshPetscCall(MatSetFromOptions(this->_mat));
- LibmeshPetscCall(MatSeqBAIJSetPreallocation(this->_mat, blocksize, n_nz/blocksize, NULL));
- LibmeshPetscCall(MatMPIBAIJSetPreallocation(this->_mat, blocksize,
- n_nz/blocksize, NULL,
- n_oz/blocksize, NULL));
}
else
#endif
@@ -196,10 +185,7 @@ void PetscMatrix::init (const numeric_index_type m_in,
// MatSetFromOptions needs to happen before Preallocation routines
// since MatSetFromOptions can change matrix type and remove incompatible
// preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
LibmeshPetscCall(MatSetFromOptions(this->_mat));
- LibmeshPetscCall(MatSeqAIJSetPreallocation(this->_mat, n_nz, NULL));
- LibmeshPetscCall(MatMPIAIJSetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
break;
case HYPRE:
@@ -209,9 +195,7 @@ void PetscMatrix::init (const numeric_index_type m_in,
// MatSetFromOptions needs to happen before Preallocation routines
// since MatSetFromOptions can change matrix type and remove incompatible
// preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
LibmeshPetscCall(MatSetFromOptions(this->_mat));
- LibmeshPetscCall(MatHYPRESetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
#else
libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
#endif
@@ -221,64 +205,74 @@ void PetscMatrix::init (const numeric_index_type m_in,
}
}
- // Make it an error for PETSc to allocate new nonzero entries during assembly
- LibmeshPetscCall(MatSetOption(this->_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
-
this->set_context ();
- this->zero ();
}
+
template
void PetscMatrix::init (const numeric_index_type m_in,
const numeric_index_type n_in,
const numeric_index_type m_l,
const numeric_index_type n_l,
- const std::vector & n_nz,
- const std::vector & n_oz,
+ const numeric_index_type nnz,
+ const numeric_index_type noz,
const numeric_index_type blocksize_in)
{
- // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
- libmesh_ignore(blocksize_in);
+ this->init_without_preallocation(m_in, n_in, m_l, n_l, blocksize_in);
- PetscInt blocksize = static_cast(blocksize_in);
+ PetscInt n_nz = static_cast(nnz);
+ PetscInt n_oz = static_cast(noz);
- // Clear initialized matrices
- if (this->initialized())
- this->clear();
+#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
+ if (blocksize > 1)
+ {
+ LibmeshPetscCall(MatSeqBAIJSetPreallocation(this->_mat, blocksize, n_nz/blocksize, NULL));
+ LibmeshPetscCall(MatMPIBAIJSetPreallocation(this->_mat, blocksize,
+ n_nz/blocksize, NULL,
+ n_oz/blocksize, NULL));
+ }
+ else
+#endif
+ {
+ switch (this->_mat_type) {
+ case AIJ:
+ LibmeshPetscCall(MatSeqAIJSetPreallocation(this->_mat, n_nz, NULL));
+ LibmeshPetscCall(MatMPIAIJSetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
+ break;
- this->_is_initialized = true;
+ case HYPRE:
+#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
+ LibmeshPetscCall(MatHYPRESetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
+#else
+ libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
+#endif
+ break;
+
+ default: libmesh_error_msg("Unsupported petsc matrix type");
+ }
+ }
+ this->finish_initialization();
+}
+
+template
+void PetscMatrix::preallocate(const numeric_index_type libmesh_dbg_var(m_l),
+ const std::vector & n_nz,
+ const std::vector & n_oz,
+ const numeric_index_type blocksize_in)
+{
// Make sure the sparsity pattern isn't empty unless the matrix is 0x0
libmesh_assert_equal_to (n_nz.size(), m_l);
libmesh_assert_equal_to (n_oz.size(), m_l);
-
- PetscInt m_global = static_cast(m_in);
- PetscInt n_global = static_cast(n_in);
- PetscInt m_local = static_cast(m_l);
- PetscInt n_local = static_cast(n_l);
-
- LibmeshPetscCall(MatCreate(this->comm().get(), &this->_mat));
- LibmeshPetscCall(MatSetSizes(this->_mat, m_local, n_local, m_global, n_global));
- LibmeshPetscCall(MatSetBlockSize(this->_mat,blocksize));
+ // Avoid unused warnings when not configured with block storage
+ libmesh_ignore(blocksize_in);
#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
+ PetscInt blocksize = static_cast(blocksize_in);
+
if (blocksize > 1)
{
- // specified blocksize, bs>1.
- // double check sizes.
- libmesh_assert_equal_to (m_local % blocksize, 0);
- libmesh_assert_equal_to (n_local % blocksize, 0);
- libmesh_assert_equal_to (m_global % blocksize, 0);
- libmesh_assert_equal_to (n_global % blocksize, 0);
-
- LibmeshPetscCall(MatSetType(this->_mat, MATBAIJ)); // Automatically chooses seqbaij or mpibaij
-
- // MatSetFromOptions needs to happen before Preallocation routines
- // since MatSetFromOptions can change matrix type and remove incompatible
- // preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
- LibmeshPetscCall(MatSetFromOptions(this->_mat));
// transform the per-entry n_nz and n_oz arrays into their block counterparts.
std::vector b_n_nz, b_n_oz;
@@ -303,13 +297,6 @@ void PetscMatrix::init (const numeric_index_type m_in,
{
switch (this->_mat_type) {
case AIJ:
- LibmeshPetscCall(MatSetType(this->_mat, MATAIJ)); // Automatically chooses seqaij or mpiaij
-
- // MatSetFromOptions needs to happen before Preallocation routines
- // since MatSetFromOptions can change matrix type and remove incompatible
- // preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
- LibmeshPetscCall(MatSetFromOptions(this->_mat));
LibmeshPetscCall(MatSeqAIJSetPreallocation (this->_mat,
0,
numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data())));
@@ -322,13 +309,6 @@ void PetscMatrix::init (const numeric_index_type m_in,
case HYPRE:
#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
- LibmeshPetscCall(MatSetType(this->_mat, MATHYPRE));
-
- // MatSetFromOptions needs to happen before Preallocation routines
- // since MatSetFromOptions can change matrix type and remove incompatible
- // preallocation
- LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
- LibmeshPetscCall(MatSetFromOptions(this->_mat));
LibmeshPetscCall(MatHYPRESetPreallocation (this->_mat,
0,
numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
@@ -343,11 +323,30 @@ void PetscMatrix::init (const numeric_index_type m_in,
}
}
+}
- // Make it an error for PETSc to allocate new nonzero entries during assembly
+template
+void PetscMatrix::finish_initialization()
+{
+ // Make it an error for PETSc to allocate new nonzero entries during assembly. For old PETSc
+ // versions this option must be set after preallocation for MPIAIJ matrices
LibmeshPetscCall(MatSetOption(this->_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
- this->set_context();
- this->zero();
+ this->_is_initialized = true;
+}
+
+template
+void PetscMatrix::init (const numeric_index_type m_in,
+ const numeric_index_type n_in,
+ const numeric_index_type m_l,
+ const numeric_index_type n_l,
+ const std::vector & n_nz,
+ const std::vector & n_oz,
+ const numeric_index_type blocksize_in)
+{
+ this->init_without_preallocation(m_in, n_in, m_l, n_l, blocksize_in);
+ this->preallocate(m_l, n_nz, n_oz, blocksize_in);
+
+ this->finish_initialization();
}
@@ -359,12 +358,18 @@ void PetscMatrix::init (const ParallelType)
const numeric_index_type m_in = this->_dof_map->n_dofs();
const numeric_index_type m_l = this->_dof_map->n_local_dofs();
- const std::vector & n_nz = this->_sp->get_n_nz();
- const std::vector & n_oz = this->_sp->get_n_oz();
+ const auto blocksize = this->_dof_map->block_size();
- PetscInt blocksize = static_cast(this->_dof_map->block_size());
+ this->init_without_preallocation(m_in, m_in, m_l, m_l, blocksize);
+ if (!this->_use_hash_table)
+ {
+ const std::vector & n_nz = this->_sp->get_n_nz();
+ const std::vector & n_oz = this->_sp->get_n_oz();
+
+ this->preallocate(m_l, n_nz, n_oz, blocksize);
+ }
- this->init(m_in, m_in, m_l, m_l, n_nz, n_oz, blocksize);
+ this->finish_initialization();
}
@@ -1267,6 +1272,45 @@ void PetscMatrix::scale(const T scale)
LibmeshPetscCall(MatScale(this->_mat, scale));
}
+template
+bool PetscMatrix::supports_hash_table() const
+{
+#if PETSC_RELEASE_LESS_THAN(3,19,0)
+ return false;
+#else
+ return true;
+#endif
+}
+
+#if PETSC_RELEASE_GREATER_EQUALS(3,23,0)
+template
+std::unique_ptr>
+PetscMatrix::copy_from_hash ()
+{
+ Mat xaij;
+ libmesh_assert(this->initialized());
+ libmesh_assert(!this->closed());
+ LibmeshPetscCall(MatDuplicate(this->_mat, MAT_DO_NOT_COPY_VALUES, &xaij));
+ LibmeshPetscCall(MatCopyHashToXAIJ(this->_mat, xaij));
+ return std::make_unique>(xaij, this->comm(), /*destroy_on_exit=*/true);
+}
+#endif
+
+template
+void
+PetscMatrix::reset_memory()
+{
+ if (this->_use_hash_table)
+#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
+ // This performs MatReset plus re-establishes the hash table
+ LibmeshPetscCall(MatResetHash(this->_mat));
+#else
+ libmesh_error_msg("Resetting hash tables not supported until PETSc version 3.23");
+#endif
+ else
+ this->reset_preallocation();
+}
+
//------------------------------------------------------------------
// Explicit instantiations
template class LIBMESH_EXPORT PetscMatrix;
diff --git a/src/numerics/petsc_matrix_base.C b/src/numerics/petsc_matrix_base.C
index 834941c99e6..df4d8acaa08 100644
--- a/src/numerics/petsc_matrix_base.C
+++ b/src/numerics/petsc_matrix_base.C
@@ -48,9 +48,10 @@ PetscMatrixBase::PetscMatrixBase(const Parallel::Communicator & comm_in) :
// for destroying it
template
PetscMatrixBase::PetscMatrixBase(Mat mat_in,
- const Parallel::Communicator & comm_in) :
+ const Parallel::Communicator & comm_in,
+ const bool destroy_on_exit) :
SparseMatrix(comm_in),
- _destroy_mat_on_exit(false)
+ _destroy_mat_on_exit(destroy_on_exit)
{
this->_mat = mat_in;
this->_is_initialized = true;
diff --git a/src/numerics/sparse_matrix.C b/src/numerics/sparse_matrix.C
index ef3f649794e..28419956f68 100644
--- a/src/numerics/sparse_matrix.C
+++ b/src/numerics/sparse_matrix.C
@@ -64,7 +64,8 @@ SparseMatrix::SparseMatrix (const Parallel::Communicator & comm_in) :
ParallelObject(comm_in),
_dof_map(nullptr),
_sp(nullptr),
- _is_initialized(false)
+ _is_initialized(false),
+ _use_hash_table(false)
{}
diff --git a/src/systems/system.C b/src/systems/system.C
index 5c4fa61021a..52754dbe0b5 100644
--- a/src/systems/system.C
+++ b/src/systems/system.C
@@ -97,7 +97,9 @@ System::System (EquationSystems & es,
_additional_data_written (false),
adjoint_already_solved (false),
_hide_output (false),
- project_with_constraints (true)
+ project_with_constraints (true),
+ _prefer_hash_table_matrix_assembly(false),
+ _require_sparsity_pattern (false)
{
}
@@ -342,18 +344,28 @@ void System::init_matrices ()
// want to attach the same matrix to the DofMap twice
if (!this->get_dof_map().is_attached(m))
this->get_dof_map().attach_matrix(m);
+
+ // If the user has already explicitly requested that this matrix use a hash table, then we
+ // always honor that
+ const bool use_hash =
+ pr.second->use_hash_table() ||
+ (this->_prefer_hash_table_matrix_assembly && pr.second->supports_hash_table());
+ pr.second->use_hash_table(use_hash);
+ if (!use_hash)
+ this->_require_sparsity_pattern = true;
}
// Compute the sparsity pattern for the current
// mesh and DOF distribution. This also updates
// additional matrices, \p DofMap now knows them
- this->get_dof_map().compute_sparsity(this->get_mesh());
+ if (this->_require_sparsity_pattern)
+ this->get_dof_map().compute_sparsity(this->get_mesh());
// Initialize matrices and set to zero
- for (auto & pr : _matrices)
+ for (auto & [name, mat] : _matrices)
{
- pr.second->init(_matrix_types[pr.first]);
- pr.second->zero();
+ mat->init(_matrix_types[name]);
+ mat->zero();
}
}
@@ -443,13 +455,16 @@ void System::reinit ()
pr.second->attach_dof_map(this->get_dof_map());
}
- // Clear the sparsity pattern
- this->get_dof_map().clear_sparsity();
+ if (this->_require_sparsity_pattern)
+ {
+ // Clear the sparsity pattern
+ this->get_dof_map().clear_sparsity();
- // Compute the sparsity pattern for the current
- // mesh and DOF distribution. This also updates
- // additional matrices, \p DofMap now knows them
- this->get_dof_map().compute_sparsity (this->get_mesh());
+ // Compute the sparsity pattern for the current
+ // mesh and DOF distribution. This also updates
+ // additional matrices, \p DofMap now knows them
+ this->get_dof_map().compute_sparsity (this->get_mesh());
+ }
// Initialize matrices and set to zero
for (auto & pr : _matrices)
diff --git a/tests/numerics/petsc_matrix_test.C b/tests/numerics/petsc_matrix_test.C
index 9835a9b56cc..9e23894ad35 100644
--- a/tests/numerics/petsc_matrix_test.C
+++ b/tests/numerics/petsc_matrix_test.C
@@ -25,6 +25,9 @@ public:
CPPUNIT_TEST(testPetscBinaryRead);
CPPUNIT_TEST(testPetscBinaryWrite);
+#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
+ CPPUNIT_TEST(testPetscCopyFromHash);
+#endif
// CPPUNIT_TEST(testPetscHDF5Read);
// CPPUNIT_TEST(testPetscHDF5Write);
@@ -99,6 +102,26 @@ public:
CPPUNIT_ASSERT_LESS(TOLERANCE, mat_to_read->l1_norm());
}
+#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
+ void testPetscCopyFromHash()
+ {
+ PetscMatrix mat(*my_comm);
+ const numeric_index_type M = my_comm->size();
+ const numeric_index_type m = 1;
+ const numeric_index_type blocksize = 1;
+ mat.use_hash_table(true);
+ mat.init_without_preallocation(M, M, m, m, blocksize);
+ mat.finish_initialization();
+
+ // We'll just write a dense row
+ for (const auto j : make_range(my_comm->size()))
+ mat.add(my_comm->rank(), j, my_comm->rank() * j);
+
+ auto copy = mat.copy_from_hash();
+ for (const auto j : make_range(my_comm->size()))
+ CPPUNIT_ASSERT_EQUAL((*copy)(my_comm->rank(), j), static_cast(my_comm->rank() * j));
+ }
+#endif
};
CPPUNIT_TEST_SUITE_REGISTRATION(PetscMatrixTest);