libMesh
Loading...
Searching...
No Matches
Public Member Functions | Static Public Member Functions | Protected Types | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions | Private Attributes | Friends | List of all members
libMesh::PetscMatrix< T > Class Template Referencefinal

This class provides a nice interface to the PETSc C-based AIJ data structures for parallel, sparse matrices. More...

#include <petsc_matrix.h>

Inheritance diagram for libMesh::PetscMatrix< T >:
[legend]

Public Member Functions

 PetscMatrix (const Parallel::Communicator &comm_in)
 Constructor; initializes the matrix to be empty, without any structure, i.e.
 
 PetscMatrix (Mat m, const Parallel::Communicator &comm_in, bool destroy_on_exit=false)
 Constructor.
 
 PetscMatrix (const Parallel::Communicator &comm_in, const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type n_nz=30, const numeric_index_type n_oz=10, const numeric_index_type blocksize=1)
 Constructor.
 
 PetscMatrix (PetscMatrix &&)=delete
 This class manages a C-style struct (Mat) manually, so we don't want to allow any automatic copy/move functions to be generated, and we can't default the destructor.
 
 PetscMatrix (const PetscMatrix &)=delete
 
PetscMatrixoperator= (PetscMatrix &&)=delete
 
virtual ~PetscMatrix ()
 
PetscMatrixoperator= (const PetscMatrix &)
 
virtual SparseMatrix< T > & operator= (const SparseMatrix< T > &v) override
 This looks like a copy assignment operator, but note that, unlike normal copy assignment operators, it is pure virtual.
 
virtual void init (const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type n_nz=30, const numeric_index_type n_oz=10, const numeric_index_type blocksize=1) override
 Initialize a PETSc matrix.
 
void init (const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const std::vector< numeric_index_type > &n_nz, const std::vector< numeric_index_type > &n_oz, const numeric_index_type blocksize=1)
 Initialize a PETSc matrix.
 
virtual void init (ParallelType=PARALLEL) override
 Initialize this matrix using the sparsity structure computed by dof_map.
 
void update_preallocation_and_zero ()
 Update the sparsity pattern based on dof_map, and set the matrix to zero.
 
void reset_preallocation ()
 Reset matrix to use the original nonzero pattern provided by users.
 
virtual void zero () override
 Set all entries to 0.
 
virtual std::unique_ptr< SparseMatrix< T > > zero_clone () const override
 
virtual std::unique_ptr< SparseMatrix< T > > clone () const override
 
virtual void zero_rows (std::vector< numeric_index_type > &rows, T diag_value=0.0) override
 Sets all row entries to 0 then puts diag_value in the diagonal entry.
 
virtual void flush () override
 For PETSc matrix , this function is similar to close but without shrinking memory.
 
virtual void set (const numeric_index_type i, const numeric_index_type j, const T value) override
 Set the element (i,j) to value.
 
virtual void add (const numeric_index_type i, const numeric_index_type j, const T value) override
 Add value to the element (i,j).
 
virtual void add_matrix (const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) override
 Add the full matrix dm to the SparseMatrix.
 
virtual void add_matrix (const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &dof_indices) override
 Same as add_matrix, but assumes the row and column maps are the same.
 
virtual void add_block_matrix (const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &brows, const std::vector< numeric_index_type > &bcols) override
 Add the full matrix dm to the SparseMatrix.
 
virtual void add_block_matrix (const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &dof_indices) override
 Same as add_block_matrix(), but assumes the row and column maps are the same.
 
virtual void add (const T a, const SparseMatrix< T > &X) override
 Compute A += a*X for scalar a, matrix X.
 
virtual void matrix_matrix_mult (SparseMatrix< T > &X, SparseMatrix< T > &Y, bool reuse=false) override
 Compute Y = A*X for matrix X.
 
virtual void add_sparse_matrix (const SparseMatrix< T > &spm, const std::map< numeric_index_type, numeric_index_type > &row_ltog, const std::map< numeric_index_type, numeric_index_type > &col_ltog, const T scalar) override
 Add scalar* spm to the rows and cols of this matrix (A): A(rows[i], cols[j]) += scalar * spm(i,j)
 
virtual T operator() (const numeric_index_type i, const numeric_index_type j) const override
 
virtual Real l1_norm () const override
 
Real frobenius_norm () const
 
virtual Real linfty_norm () const override
 
virtual void print_personal (std::ostream &os=libMesh::out) const override
 Print the contents of the matrix to the screen with the PETSc viewer.
 
virtual void print_matlab (const std::string &name="") const override
 Print the contents of the matrix in Matlab's sparse matrix format.
 
virtual void print_petsc_binary (const std::string &filename) const override
 Write the contents of the matrix to a file in PETSc's binary sparse matrix format.
 
virtual void print_petsc_hdf5 (const std::string &filename) const override
 Write the contents of the matrix to a file in PETSc's HDF5 sparse matrix format.
 
virtual void read_petsc_binary (const std::string &filename) override
 Read the contents of the matrix from a file in PETSc's binary sparse matrix format.
 
virtual void read_petsc_hdf5 (const std::string &filename) override
 Read the contents of the matrix from a file in PETSc's HDF5 sparse matrix format.
 
virtual void get_diagonal (NumericVector< T > &dest) const override
 Copies the diagonal part of the matrix into dest.
 
virtual void get_transpose (SparseMatrix< T > &dest) const override
 Copies the transpose of the matrix into dest, which may be *this.
 
virtual void get_row (numeric_index_type i, std::vector< numeric_index_type > &indices, std::vector< T > &values) const override
 Get a row from the matrix.
 
virtual void create_submatrix_nosort (SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) const override
 Similar to the create_submatrix function, this function creates a submatrix which is defined by the indices given in the rows and cols vectors.
 
virtual void scale (const T scale) override
 Scales all elements of this matrix by scale.
 
std::unique_ptr< PetscMatrix< T > > copy_from_hash ()
 Creates a copy of the current hash table matrix and then performs assembly.
 
virtual bool supports_hash_table () const override
 
virtual void restore_original_nonzero_pattern () override
 Reset the memory storage of the matrix.
 
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)
 Perform matrix initialization steps sans preallocation.
 
void finish_initialization ()
 Finish up the initialization process.
 
virtual SolverPackage solver_package () override
 
Mat mat ()
 
Mat mat () const
 
virtual void clear () noexcept override
 clear() is called from the destructor, so it should not throw.
 
void set_destroy_mat_on_exit (bool destroy=true)
 If set to false, we don't delete the Mat on destruction and allow instead for PETSc to manage it.
 
void swap (PetscMatrixBase< T > &)
 Swaps the internal data pointers of two PetscMatrices, no actual values are swapped.
 
void set_context ()
 Set the context (ourself) for _mat.
 
virtual numeric_index_type m () const override
 
virtual numeric_index_type local_m () const final
 Get the number of rows owned by this process.
 
virtual numeric_index_type n () const override
 
virtual numeric_index_type local_n () const final
 Get the number of columns owned by this process.
 
virtual numeric_index_type row_start () const override
 
virtual numeric_index_type row_stop () const override
 
virtual numeric_index_type col_start () const override
 
virtual numeric_index_type col_stop () const override
 
virtual void close () override
 Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across processors.
 
virtual bool closed () const override
 
virtual bool initialized () const
 
void attach_dof_map (const DofMap &dof_map)
 Set a pointer to the DofMap to use.
 
void attach_sparsity_pattern (const SparsityPattern::Build &sp)
 Set a pointer to a sparsity pattern to use.
 
virtual bool need_full_sparsity_pattern () const
 
virtual bool require_sparsity_pattern () const
 
virtual void update_sparsity_pattern (const SparsityPattern::Graph &)
 Updates the matrix sparsity pattern.
 
Real l1_norm_diff (const SparseMatrix< T > &other_mat) const
 
virtual std::size_t n_nonzeros () const
 
void print (std::ostream &os=libMesh::out, const bool sparse=false) const
 Print the contents of the matrix to the screen in a uniform style, regardless of matrix/solver package being used.
 
void print (const std::string &filename) const
 Print the contents of the matrix to a file, with a file format depending on the extension of filename.
 
void print (std::ostream &os, const bool sparse) const
 
virtual void print_coreform_hdf5 (const std::string &filename, const std::string &groupname="extraction") const
 Print the contents of the matrix to a file, with the HDF5 sparse matrix format used by CoreForm, putting CSR sparse matrix data in the group given by groupname.
 
virtual void read (const std::string &filename)
 Read the contents of the matrix from a file, with the file format inferred from the extension of filename.
 
virtual void read_coreform_hdf5 (const std::string &filename, const std::string &groupname="extraction")
 Read the contents of the matrix from a file, with the HDF5 sparse matrix format used by CoreForm, expecing sparse matrix data in the group given by groupname.
 
virtual void read_matlab (const std::string &filename)
 Read the contents of the matrix from the Matlab-script sparse matrix format used by PETSc.
 
virtual void create_submatrix (SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) const
 This function creates a matrix called "submatrix" which is defined by the row and column indices given in the "rows" and "cols" entries.
 
virtual void reinit_submatrix (SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) const
 This function is similar to the one above, but it allows you to reuse the existing sparsity pattern of "submatrix" instead of reallocating it again.
 
void vector_mult (NumericVector< T > &dest, const NumericVector< T > &arg) const
 Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest.
 
void vector_mult_add (NumericVector< T > &dest, const NumericVector< T > &arg) const
 Multiplies the matrix by the NumericVector arg and adds the result to the NumericVector dest.
 
void use_hash_table (bool use_hash)
 Sets whether to use hash table assembly.
 
bool use_hash_table () const
 
const Parallel::Communicatorcomm () const
 
processor_id_type n_processors () const
 
processor_id_type processor_id () const
 

Static Public Member Functions

static PetscMatrixBase< T > * get_context (Mat mat, const TIMPI::Communicator &comm)
 
static std::unique_ptr< SparseMatrix< T > > build (const Parallel::Communicator &comm, const SolverPackage solver_package=libMesh::default_solver_package(), const MatrixBuildType matrix_build_type=MatrixBuildType::AUTOMATIC)
 Builds a SparseMatrix<T> using the linear solver package specified by solver_package.
 
static std::string get_info ()
 Gets a string containing the reference information.
 
static void print_info (std::ostream &out_stream=libMesh::out)
 Prints the reference information, by default to libMesh::out.
 
static unsigned int n_objects ()
 Prints the number of outstanding (created, but not yet destroyed) objects.
 
static void enable_print_counter_info ()
 Methods to enable/disable the reference counter output from print_info().
 
static void disable_print_counter_info ()
 

Protected Types

typedef std::map< std::string, std::pair< unsigned int, unsigned int > > Counts
 Data structure to log the information.
 

Protected Member Functions

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)
 
virtual void _get_submatrix (SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols, const bool reuse_submatrix) const override
 This function either creates or re-initializes a matrix called submatrix which is defined by the indices given in the rows and cols vectors.
 
void _petsc_viewer (const std::string &filename, PetscViewerType viewertype, PetscFileMode filemode)
 
void increment_constructor_count (const std::string &name) noexcept
 Increments the construction counter.
 
void increment_destructor_count (const std::string &name) noexcept
 Increments the destruction counter.
 

Protected Attributes

PetscMatrixType _mat_type
 
Mat _mat
 PETSc matrix datatype to store values.
 
bool _destroy_mat_on_exit
 This boolean value should only be set to false for the constructor which takes a PETSc Mat object.
 
DofMap const * _dof_map
 The DofMap object associated with this object.
 
SparsityPattern::Build const * _sp
 The sparsity pattern associated with this object.
 
bool _is_initialized
 Flag indicating whether or not the matrix has been initialized.
 
bool _use_hash_table
 Flag indicating whether the matrix is assembled using a hash table.
 
const Parallel::Communicator_communicator
 

Static Protected Attributes

static Counts _counts
 Actually holds the data.
 
static Threads::atomic< unsigned int_n_objects
 The number of objects.
 
static Threads::spin_mutex _mutex
 Mutual exclusion object to enable thread-safe reference counting.
 
static bool _enable_print_counter = true
 Flag to control whether reference count information is printed when print_info is called.
 

Private Member Functions

template<NormType N>
Real norm () const
 

Private Attributes

std::mutex _petsc_matrix_mutex
 
Threads::spin_mutex _petsc_matrix_mutex
 

Friends

class ::PetscMatrixTest
 

Detailed Description

template<typename T>
class libMesh::PetscMatrix< T >

This class provides a nice interface to the PETSc C-based AIJ data structures for parallel, sparse matrices.

All overridden virtual functions are documented in sparse_matrix.h.

Author
Benjamin S. Kirk
Date
2002

SparseMatrix interface to PETSc Mat.

Definition at line 61 of file petsc_matrix.h.

Member Typedef Documentation

◆ Counts

typedef std::map<std::string, std::pair<unsigned int, unsigned int> > libMesh::ReferenceCounter::Counts
protectedinherited

Data structure to log the information.

The log is identified by the class name.

Definition at line 119 of file reference_counter.h.

Constructor & Destructor Documentation

◆ PetscMatrix() [1/5]

template<typename T >
libMesh::PetscMatrix< T >::PetscMatrix ( const Parallel::Communicator comm_in)
explicit

Constructor; initializes the matrix to be empty, without any structure, i.e.

the matrix is not usable at all. This constructor is therefore only useful for matrices which are members of a class. All other matrices should be created at a point in the data flow where all necessary information is available.

You have to initialize the matrix before usage with init(...).

Definition at line 86 of file petsc_matrix.C.

86 :
88{
89}
This class provides a nice interface to the PETSc C-based data structures for parallel,...
PetscMatrixType _mat_type

◆ PetscMatrix() [2/5]

template<typename T >
libMesh::PetscMatrix< T >::PetscMatrix ( Mat  m,
const Parallel::Communicator comm_in,
bool  destroy_on_exit = false 
)
explicit

Constructor.

Creates a PetscMatrix assuming you already have a 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.

Definition at line 96 of file petsc_matrix.C.

98 :
99 PetscMatrixBase<T>(mat_in, comm_in, destroy_on_exit)
100{
101 MatType mat_type;
102 LibmeshPetscCall(MatGetType(mat_in, &mat_type));
103 PetscBool is_hypre;
104 LibmeshPetscCall(PetscStrcmp(mat_type, MATHYPRE, &is_hypre));
105 if (is_hypre == PETSC_TRUE)
107 else
108 _mat_type = AIJ;
109}

References libMesh::PetscMatrix< T >::_mat_type, libMesh::AIJ, and libMesh::HYPRE.

◆ PetscMatrix() [3/5]

template<typename T >
libMesh::PetscMatrix< T >::PetscMatrix ( const Parallel::Communicator comm_in,
const numeric_index_type  m,
const numeric_index_type  n,
const numeric_index_type  m_l,
const numeric_index_type  n_l,
const numeric_index_type  n_nz = 30,
const numeric_index_type  n_oz = 10,
const numeric_index_type  blocksize = 1 
)
explicit

Constructor.

Creates and initializes a PetscMatrix with the given structure. See init(...) for a description of the parameters.

Definition at line 115 of file petsc_matrix.C.

122 :
124{
125 this->init(m_in, n_in, m_l, n_l, n_nz, n_oz, blocksize_in);
126}
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type n_nz=30, const numeric_index_type n_oz=10, const numeric_index_type blocksize=1) override
Initialize a PETSc matrix.

References libMesh::PetscMatrix< T >::init().

◆ PetscMatrix() [4/5]

template<typename T >
libMesh::PetscMatrix< T >::PetscMatrix ( PetscMatrix< T > &&  )
delete

This class manages a C-style struct (Mat) manually, so we don't want to allow any automatic copy/move functions to be generated, and we can't default the destructor.

◆ PetscMatrix() [5/5]

template<typename T >
libMesh::PetscMatrix< T >::PetscMatrix ( const PetscMatrix< T > &  )
delete

◆ ~PetscMatrix()

template<typename T >
libMesh::PetscMatrix< T >::~PetscMatrix ( )
virtualdefault

Member Function Documentation

◆ _get_submatrix()

template<typename T >
void libMesh::PetscMatrix< T >::_get_submatrix ( SparseMatrix< T > &  submatrix,
const std::vector< numeric_index_type > &  rows,
const std::vector< numeric_index_type > &  cols,
const bool  reuse_submatrix 
) const
overrideprotectedvirtual

This function either creates or re-initializes a matrix called submatrix which is defined by the indices given in the rows and cols vectors.

This function is implemented in terms of MatGetSubMatrix(). The reuse_submatrix parameter determines whether or not PETSc will treat submatrix as one which has already been used (had memory allocated) or as a new matrix.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 802 of file petsc_matrix.C.

806{
807 if (!this->closed())
808 {
809 libmesh_deprecated();
810 libmesh_warning("The matrix must be assembled before calling PetscMatrix::create_submatrix().\n"
811 "Please update your code, as this warning will become an error in a future release.");
812 const_cast<PetscMatrix<T> *>(this)->close();
813 }
814
815 semiparallel_only();
816
817 // Make sure the SparseMatrix passed in is really a PetscMatrix
818 PetscMatrix<T> * petsc_submatrix = cast_ptr<PetscMatrix<T> *>(&submatrix);
819
820 // If we're not reusing submatrix and submatrix is already initialized
821 // then we need to clear it, otherwise we get a memory leak.
822 if (!reuse_submatrix && submatrix.initialized())
823 submatrix.clear();
824
825 // Construct row and column index sets.
826 WrappedPetsc<IS> isrow;
827 LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
828 cast_int<PetscInt>(rows.size()),
829 numeric_petsc_cast(rows.data()),
830 PETSC_USE_POINTER,
831 isrow.get()));
832
833 WrappedPetsc<IS> iscol;
834 LibmeshPetscCall(ISCreateGeneral(this->comm().get(),
835 cast_int<PetscInt>(cols.size()),
836 numeric_petsc_cast(cols.data()),
837 PETSC_USE_POINTER,
838 iscol.get()));
839
840 // Extract submatrix
841 LibmeshPetscCall(LibMeshCreateSubMatrix(this->_mat,
842 isrow,
843 iscol,
844 (reuse_submatrix ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX),
845 (&petsc_submatrix->_mat)));
846
847 // Specify that the new submatrix is initialized and close it.
848 petsc_submatrix->_is_initialized = true;
849 petsc_submatrix->close();
850}
const Parallel::Communicator & comm() const
virtual void close() override
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
virtual bool closed() const override
Mat _mat
PETSc matrix datatype to store values.
This class provides a nice interface to the PETSc C-based AIJ data structures for parallel,...
virtual bool initialized() const
virtual void clear()=0
Restores the SparseMatrix<T> to a pristine state.
bool _is_initialized
Flag indicating whether or not the matrix has been initialized.
const Elem & get(const ElemType type_in)
PetscInt * numeric_petsc_cast(const numeric_index_type *p)

References libMesh::SparseMatrix< T >::_is_initialized, libMesh::PetscMatrixBase< T >::_mat, libMesh::SparseMatrix< T >::clear(), libMesh::PetscMatrixBase< T >::close(), libMesh::closed(), libMesh::WrappedPetsc< T >::get(), libMesh::SparseMatrix< T >::initialized(), and libMesh::numeric_petsc_cast().

◆ _petsc_viewer()

template<typename T >
void libMesh::PetscMatrix< T >::_petsc_viewer ( const std::string &  filename,
PetscViewerType  viewertype,
PetscFileMode  filemode 
)
protected

Definition at line 661 of file petsc_matrix.C.

664{
665 parallel_object_only();
666
667 // We'll get matrix sizes from the file, but we need to at least
668 // have a Mat object
669 if (!this->initialized())
670 {
671 LibmeshPetscCall(MatCreate(this->comm().get(), &this->_mat));
672 this->_is_initialized = true;
673 }
674
675 PetscViewer viewer;
676 LibmeshPetscCall(PetscViewerCreate(this->comm().get(), &viewer));
677 LibmeshPetscCall(PetscViewerSetType(viewer, viewertype));
678 LibmeshPetscCall(PetscViewerSetFromOptions(viewer));
679 LibmeshPetscCall(PetscViewerFileSetMode(viewer, filemode));
680 LibmeshPetscCall(PetscViewerFileSetName(viewer, filename.c_str()));
681 if (filemode == FILE_MODE_READ)
682 LibmeshPetscCall(MatLoad(this->_mat, viewer));
683 else
684 LibmeshPetscCall(MatView(this->_mat, viewer));
685 LibmeshPetscCall(PetscViewerDestroy(&viewer));
686}
bool initialized()
Checks that library initialization has been done.
Definition libmesh.C:324

References libMesh::initialized().

Referenced by libMesh::PetscMatrix< T >::print_petsc_binary(), and libMesh::PetscMatrix< T >::print_petsc_hdf5().

◆ add() [1/2]

template<typename T >
void libMesh::PetscMatrix< T >::add ( const numeric_index_type  i,
const numeric_index_type  j,
const T  value 
)
overridevirtual

Add value to the element (i,j).

Throws an error if the entry does not exist. Zero values can be "added" to non-existent entries.

Implements libMesh::SparseMatrix< T >.

Definition at line 989 of file petsc_matrix.C.

992{
993 libmesh_assert (this->initialized());
994
995 PetscInt i_val=i, j_val=j;
996
997 PetscScalar petsc_value = static_cast<PetscScalar>(value);
998 std::scoped_lock lock(this->_petsc_matrix_mutex);
999 LibmeshPetscCall(MatSetValues(this->_mat, 1, &i_val, 1, &j_val,
1000 &petsc_value, ADD_VALUES));
1001}
std::mutex _petsc_matrix_mutex
libmesh_assert(ctx)
static const bool value
Definition xdr_io.C:55

References libMesh::initialized(), libMesh::libmesh_assert(), and value.

Referenced by PetscMatrixTest::testPetscCopyFromHash().

◆ add() [2/2]

template<typename T >
void libMesh::PetscMatrix< T >::add ( const T  a,
const SparseMatrix< T > &  X 
)
overridevirtual

Compute A += a*X for scalar a, matrix X.

Note
The matrices A and X need to have the same nonzero pattern, otherwise PETSc will crash!
It is advisable to not only allocate appropriate memory with init(), but also explicitly zero the terms of this whenever you add a non-zero value to X.
X will be closed, if not already done, before performing any work.

Implements libMesh::SparseMatrix< T >.

Definition at line 1019 of file petsc_matrix.C.

1020{
1021 libmesh_assert (this->initialized());
1022
1023 // sanity check. but this cannot avoid
1024 // crash due to incompatible sparsity structure...
1025 libmesh_assert_equal_to (this->m(), X_in.m());
1026 libmesh_assert_equal_to (this->n(), X_in.n());
1027
1028 PetscScalar a = static_cast<PetscScalar> (a_in);
1029 const PetscMatrix<T> * X = cast_ptr<const PetscMatrix<T> *> (&X_in);
1030
1031 libmesh_assert (X);
1032
1033 // the matrix from which we copy the values has to be assembled/closed
1034 libmesh_assert(X->closed());
1035
1036 semiparallel_only();
1037
1038 LibmeshPetscCall(MatAXPY(this->_mat, a, X->_mat, DIFFERENT_NONZERO_PATTERN));
1039}
virtual numeric_index_type n() const override
virtual numeric_index_type m() const override

References libMesh::PetscMatrixBase< T >::_mat, libMesh::PetscMatrixBase< T >::closed(), libMesh::initialized(), libMesh::libmesh_assert(), libMesh::SparseMatrix< T >::m(), and libMesh::SparseMatrix< T >::n().

◆ add_block_matrix() [1/2]

template<typename T >
void libMesh::PetscMatrix< T >::add_block_matrix ( const DenseMatrix< T > &  dm,
const std::vector< numeric_index_type > &  brows,
const std::vector< numeric_index_type > &  bcols 
)
overridevirtual

Add the full matrix dm to the SparseMatrix.

This is useful for adding an element matrix at assembly time. The matrix is assumed blocked, and brow, bcol correspond to the block row and column indices.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 761 of file petsc_matrix.C.

764{
765 libmesh_assert (this->initialized());
766
767 const numeric_index_type n_brows =
768 cast_int<numeric_index_type>(brows.size());
769 const numeric_index_type n_bcols =
770 cast_int<numeric_index_type>(bcols.size());
771
772#ifndef NDEBUG
773 const numeric_index_type n_rows =
774 cast_int<numeric_index_type>(dm.m());
775 const numeric_index_type n_cols =
776 cast_int<numeric_index_type>(dm.n());
777 const numeric_index_type blocksize = n_rows / n_brows;
778
779 libmesh_assert_equal_to (n_cols / n_bcols, blocksize);
780 libmesh_assert_equal_to (blocksize*n_brows, n_rows);
781 libmesh_assert_equal_to (blocksize*n_bcols, n_cols);
782
783 PetscInt petsc_blocksize;
784 LibmeshPetscCall(MatGetBlockSize(this->_mat, &petsc_blocksize));
785 libmesh_assert_equal_to (blocksize, static_cast<numeric_index_type>(petsc_blocksize));
786#endif
787
788 std::scoped_lock lock(this->_petsc_matrix_mutex);
789 // These casts are required for PETSc <= 2.1.5
790 LibmeshPetscCall(MatSetValuesBlocked(this->_mat,
791 n_brows, numeric_petsc_cast(brows.data()),
792 n_bcols, numeric_petsc_cast(bcols.data()),
793 pPS(const_cast<T*>(dm.get_values().data())),
794 ADD_VALUES));
795}
std::vector< T > & get_values()
PetscScalar * pPS(T *ptr)
dof_id_type numeric_index_type
Definition id_types.h:99

References libMesh::DenseMatrix< T >::get_values(), libMesh::initialized(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::numeric_petsc_cast(), and libMesh::pPS().

Referenced by libMesh::PetscMatrix< T >::add_block_matrix().

◆ add_block_matrix() [2/2]

template<typename T >
virtual void libMesh::PetscMatrix< T >::add_block_matrix ( const DenseMatrix< T > &  dm,
const std::vector< numeric_index_type > &  dof_indices 
)
inlineoverridevirtual

Same as add_block_matrix(), but assumes the row and column maps are the same.

Thus the matrix dm must be square.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 197 of file petsc_matrix.h.

199 { this->add_block_matrix (dm, dof_indices, dof_indices); }
virtual void add_block_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &brows, const std::vector< numeric_index_type > &bcols) override
Add the full matrix dm to the SparseMatrix.

References libMesh::PetscMatrix< T >::add_block_matrix().

◆ add_matrix() [1/2]

template<typename T >
void libMesh::PetscMatrix< T >::add_matrix ( const DenseMatrix< T > &  dm,
const std::vector< numeric_index_type > &  dof_indices 
)
overridevirtual

Same as add_matrix, but assumes the row and column maps are the same.

Thus the matrix dm must be square.

Implements libMesh::SparseMatrix< T >.

Definition at line 1006 of file petsc_matrix.C.

1008{
1009 this->add_matrix (dm, dof_indices, dof_indices);
1010}
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) override
Add the full matrix dm to the SparseMatrix.

◆ add_matrix() [2/2]

template<typename T >
void libMesh::PetscMatrix< T >::add_matrix ( const DenseMatrix< T > &  dm,
const std::vector< numeric_index_type > &  rows,
const std::vector< numeric_index_type > &  cols 
)
overridevirtual

Add the full matrix dm to the SparseMatrix.

This is useful for adding an element matrix at assembly time.

Implements libMesh::SparseMatrix< T >.

Definition at line 735 of file petsc_matrix.C.

738{
739 libmesh_assert (this->initialized());
740
741 const numeric_index_type n_rows = dm.m();
742 const numeric_index_type n_cols = dm.n();
743
744 libmesh_assert_equal_to (rows.size(), n_rows);
745 libmesh_assert_equal_to (cols.size(), n_cols);
746
747 std::scoped_lock lock(this->_petsc_matrix_mutex);
748 LibmeshPetscCall(MatSetValues(this->_mat,
749 n_rows, numeric_petsc_cast(rows.data()),
750 n_cols, numeric_petsc_cast(cols.data()),
751 pPS(const_cast<T*>(dm.get_values().data())),
752 ADD_VALUES));
753}

References libMesh::DenseMatrix< T >::get_values(), libMesh::initialized(), libMesh::libmesh_assert(), libMesh::DenseMatrixBase< T >::m(), libMesh::DenseMatrixBase< T >::n(), libMesh::numeric_petsc_cast(), and libMesh::pPS().

◆ add_sparse_matrix()

template<typename T >
void libMesh::PetscMatrix< T >::add_sparse_matrix ( const SparseMatrix< T > &  spm,
const std::map< numeric_index_type, numeric_index_type > &  row_ltog,
const std::map< numeric_index_type, numeric_index_type > &  col_ltog,
const T  scalar 
)
overridevirtual

Add scalar* spm to the rows and cols of this matrix (A): A(rows[i], cols[j]) += scalar * spm(i,j)

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1073 of file petsc_matrix.C.

1077{
1078 // size of spm is usually greater than row_ltog and col_ltog in parallel as the indices are owned by the processor
1079 // also, we should allow adding certain parts of spm to _mat
1080 libmesh_assert_greater_equal(spm.m(), row_ltog.size());
1081 libmesh_assert_greater_equal(spm.n(), col_ltog.size());
1082
1083 // make sure matrix has larger size than spm
1084 libmesh_assert_greater_equal(this->m(), spm.m());
1085 libmesh_assert_greater_equal(this->n(), spm.n());
1086
1087 if (!this->closed())
1088 this->close();
1089
1090 auto pscm = cast_ptr<const PetscMatrix<T> *>(&spm);
1091
1092 PetscInt ncols = 0;
1093
1094 const PetscInt * lcols;
1095 const PetscScalar * vals;
1096
1097 std::vector<PetscInt> gcols;
1098 std::vector<PetscScalar> values;
1099
1100 for (auto ltog : row_ltog)
1101 {
1102 PetscInt grow[] = {static_cast<PetscInt>(ltog.second)}; // global row index
1103
1104 LibmeshPetscCall(MatGetRow(pscm->_mat, static_cast<PetscInt>(ltog.first), &ncols, &lcols, &vals));
1105
1106 // get global indices (gcols) from lcols, increment values = vals*scalar
1107 gcols.resize(ncols);
1108 values.resize(ncols);
1109 for (auto i : index_range(gcols))
1110 {
1111 gcols[i] = libmesh_map_find(col_ltog, lcols[i]);
1112 values[i] = PS(scalar) * vals[i];
1113 }
1114
1115 LibmeshPetscCall(MatSetValues(this->_mat, 1, grow, ncols, gcols.data(), values.data(), ADD_VALUES));
1116 LibmeshPetscCall(MatRestoreRow(pscm->_mat, static_cast<PetscInt>(ltog.first), &ncols, &lcols, &vals));
1117 }
1118 // Note: We are not closing the matrix because it is expensive to do so when adding multiple sparse matrices.
1119 // Remember to manually close the matrix once all changes to the matrix have been made.
1120}
virtual numeric_index_type n() const =0
virtual numeric_index_type m() const =0
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition int_range.h:153
PetscScalar PS(T val)

References libMesh::closed(), libMesh::index_range(), libMesh::SparseMatrix< T >::m(), libMesh::SparseMatrix< T >::n(), and libMesh::PS().

◆ attach_dof_map()

template<typename T >
void libMesh::SparseMatrix< T >::attach_dof_map ( const DofMap dof_map)
inherited

Set a pointer to the DofMap to use.

If a separate sparsity pattern is not being used, use the one from the DofMap.

The lifetime of dof_map must exceed the lifetime of this.

Definition at line 100 of file sparse_matrix.C.

101{
102 _dof_map = &dof_map;
103 if (!_sp)
104 _sp = dof_map.get_sparsity_pattern();
105}
DofMap const * _dof_map
The DofMap object associated with this object.
SparsityPattern::Build const * _sp
The sparsity pattern associated with this object.
const SparsityPattern::Graph & get_sparsity_pattern() const
Rows of sparse matrix indices, indexed by the offset from the first DoF on this processor.

References libMesh::DofMap::get_sparsity_pattern().

Referenced by libMesh::Problem_Interface::computeJacobian(), DMlibMeshJacobian(), and libMesh::DofMap::update_sparsity_pattern().

◆ attach_sparsity_pattern()

template<typename T >
void libMesh::SparseMatrix< T >::attach_sparsity_pattern ( const SparsityPattern::Build sp)
inherited

Set a pointer to a sparsity pattern to use.

Useful in cases where a matrix requires a wider (or for efficiency narrower) pattern than most matrices in the system, or in cases where no system sparsity pattern is being calculated by the DofMap.

The lifetime of sp must exceed the lifetime of this.

Definition at line 110 of file sparse_matrix.C.

111{
112 _sp = &sp;
113}

Referenced by libMesh::DofMap::update_sparsity_pattern().

◆ build()

template<typename T >
std::unique_ptr< SparseMatrix< T > > libMesh::SparseMatrix< T >::build ( const Parallel::Communicator comm,
const SolverPackage  solver_package = libMesh::default_solver_package(),
const MatrixBuildType  matrix_build_type = MatrixBuildType::AUTOMATIC 
)
staticinherited

Builds a SparseMatrix<T> using the linear solver package specified by solver_package.

Definition at line 193 of file sparse_matrix.C.

196{
197 // Avoid unused parameter warnings when no solver packages are enabled.
199
200 if (matrix_build_type == MatrixBuildType::DIAGONAL)
201 return std::make_unique<DiagonalMatrix<T>>(comm);
202
203 // Build the appropriate vector
204 switch (solver_package)
205 {
206
207#ifdef LIBMESH_HAVE_LASPACK
208 case LASPACK_SOLVERS:
209 return std::make_unique<LaspackMatrix<T>>(comm);
210#endif
211
212
213#ifdef LIBMESH_HAVE_PETSC
214 case PETSC_SOLVERS:
215 return std::make_unique<PetscMatrix<T>>(comm);
216#endif
217
218
219#ifdef LIBMESH_TRILINOS_HAVE_EPETRA
220 case TRILINOS_SOLVERS:
221 return std::make_unique<EpetraMatrix<T>>(comm);
222#endif
223
224
225#ifdef LIBMESH_HAVE_EIGEN
226 case EIGEN_SOLVERS:
227 return std::make_unique<EigenSparseMatrix<T>>(comm);
228#endif
229
230 default:
231 libmesh_error_msg("ERROR: Unrecognized solver package: " << solver_package);
232 }
233}
virtual SolverPackage solver_package()=0
void libmesh_ignore(const Args &...)

References libMesh::DIAGONAL, libMesh::EIGEN_SOLVERS, libMesh::LASPACK_SOLVERS, libMesh::libmesh_ignore(), libMesh::PETSC_SOLVERS, and libMesh::TRILINOS_SOLVERS.

Referenced by libMesh::CondensedEigenSystem::add_matrices(), libMesh::System::add_matrix(), libMesh::RBConstruction::allocate_data_structures(), libMesh::TransientRBConstruction::allocate_data_structures(), libMesh::CondensedEigenSystem::copy_super_to_sub(), libMesh::StaticCondensation::init(), main(), libMesh::DofMap::process_mesh_constraint_rows(), ConstraintOperatorTest::test1DCoarseningNewNodes(), ConstraintOperatorTest::test1DCoarseningOperator(), ConstraintOperatorTest::testCoreform(), ConnectedComponentsTest::testEdge(), SystemsTest::testProjectMatrix1D(), SystemsTest::testProjectMatrix2D(), and SystemsTest::testProjectMatrix3D().

◆ clear()

template<typename T >
void libMesh::PetscMatrixBase< T >::clear ( )
overridevirtualnoexceptinherited

clear() is called from the destructor, so it should not throw.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 74 of file petsc_matrix_base.C.

75{
76 if ((this->initialized()) && (this->_destroy_mat_on_exit))
77 {
78 exceptionless_semiparallel_only();
79
80 // If we encounter an error here, print a warning but otherwise
81 // keep going since we may be recovering from an exception.
82 PetscErrorCode ierr = MatDestroy (&_mat);
83 if (ierr)
84 libmesh_warning("Warning: MatDestroy returned a non-zero error code which we ignored.");
85
86 this->_is_initialized = false;
87 }
88}
bool _destroy_mat_on_exit
This boolean value should only be set to false for the constructor which takes a PETSc Mat object.

References libMesh::initialized().

Referenced by libMesh::StaticCondensation::clear().

◆ clone()

template<typename T >
std::unique_ptr< SparseMatrix< T > > libMesh::PetscMatrix< T >::clone ( ) const
overridevirtual
Returns
A smart pointer to a copy of this matrix.
Note
This must be overridden in the derived classes.

Implements libMesh::SparseMatrix< T >.

Definition at line 454 of file petsc_matrix.C.

455{
456 libmesh_error_msg_if(!this->closed(), "Matrix must be closed before it can be cloned!");
457
458 // Copy the nonzero pattern and numerical values
459 Mat copy;
460 LibmeshPetscCall(MatDuplicate(this->_mat, MAT_COPY_VALUES, &copy));
461
462 // Call wrapping PetscMatrix constructor, have it take over
463 // ownership.
464 auto ret = std::make_unique<PetscMatrix<T>>(copy, this->comm());
465 ret->set_destroy_mat_on_exit(true);
466
467 return ret;
468}

References libMesh::closed().

◆ close()

template<typename T >
void libMesh::PetscMatrixBase< T >::close ( )
overridevirtualinherited

Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across processors.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 226 of file petsc_matrix_base.C.

227{
228 semiparallel_only();
229
230 // BSK - 1/19/2004
231 // strictly this check should be OK, but it seems to
232 // fail on matrix-free matrices. Do they falsely
233 // state they are assembled? Check with the developers...
234 // if (this->closed())
235 // return;
236
237 MatAssemblyBeginEnd(this->comm(), this->_mat, MAT_FINAL_ASSEMBLY);
238}

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::PetscMatrix< T >::_get_submatrix(), libMesh::PetscLinearSolver< Number >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscMatrix< T >::create_submatrix_nosort(), DMlibMeshJacobian(), libMesh::PetscMatrix< T >::get_transpose(), libMesh::PetscPreconditioner< T >::set_hypre_ads_data(), and libMesh::PetscPreconditioner< T >::set_hypre_ams_data().

◆ closed()

template<typename T >
bool libMesh::PetscMatrixBase< T >::closed ( ) const
overridevirtualinherited
Returns
true if the matrix has been assembled.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 241 of file petsc_matrix_base.C.

242{
243 libmesh_assert (this->initialized());
244
245 PetscBool assembled;
246
247 LibmeshPetscCall(MatAssembled(this->_mat, &assembled));
248
249 return (assembled == PETSC_TRUE);
250}

References libMesh::initialized(), and libMesh::libmesh_assert().

Referenced by libMesh::PetscMatrix< T >::add(), and libMesh::PetscMatrix< T >::matrix_matrix_mult().

◆ col_start()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::col_start ( ) const
overridevirtualinherited
Returns
The index of the first matrix column owned by this processor.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 202 of file petsc_matrix_base.C.

203{
204 libmesh_assert (this->initialized());
205
206 PetscInt start=0, stop=0;
207
208 LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
209
210 return static_cast<numeric_index_type>(start);
211}
void stop(const char *file, int line, const char *date, const char *time)

References libMesh::initialized(), and libMesh::libmesh_assert().

◆ col_stop()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::col_stop ( ) const
overridevirtualinherited
Returns
The index of the last matrix column (+1) owned by this processor.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 214 of file petsc_matrix_base.C.

215{
216 libmesh_assert (this->initialized());
217
218 PetscInt start=0, stop=0;
219
220 LibmeshPetscCall(MatGetOwnershipRangeColumn(this->_mat, &start, &stop));
221
222 return static_cast<numeric_index_type>(stop);
223}

References libMesh::initialized(), and libMesh::libmesh_assert().

◆ comm()

const Parallel::Communicator & libMesh::ParallelObject::comm ( ) const
inlineinherited
Returns
A reference to the Parallel::Communicator object used by this mesh.

Definition at line 97 of file parallel_object.h.

98 { return _communicator; }
const Parallel::Communicator & _communicator

References libMesh::ParallelObject::_communicator.

Referenced by libMesh::__libmesh_petsc_diff_solver_jacobian(), libMesh::__libmesh_petsc_diff_solver_monitor(), libMesh::__libmesh_petsc_diff_solver_residual(), libMesh::ExactSolution::_compute_error(), libMesh::UniformRefinementEstimator::_estimate_error(), libMesh::Partitioner::_find_global_index_by_pid_map(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_get_diagonal(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult(), libMesh::SlepcEigenSolver< T >::_petsc_shell_matrix_mult(), libMesh::PetscLinearSolver< T >::_petsc_shell_matrix_mult_add(), libMesh::DofMap::add_constraints_to_send_list(), add_cube_convex_hull_to_mesh(), libMesh::PetscDMWrapper::add_dofs_helper(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::TransientRBConstruction::add_IC_to_RB_space(), libMesh::RBEIMEvaluation::add_interpolation_data(), libMesh::CondensedEigenSystem::add_matrices(), libMesh::EigenSystem::add_matrices(), libMesh::System::add_matrix(), libMesh::System::add_matrix(), libMesh::System::add_matrix(), libMesh::RBConstruction::add_scaled_matrix_and_vector(), libMesh::System::add_vector(), libMesh::MeshTools::Modification::all_tri(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::DofMap::allgather_recursive_constraints(), libMesh::RBConstruction::allocate_data_structures(), libMesh::TransientRBConstruction::allocate_data_structures(), libMesh::TransientRBConstruction::assemble_affine_expansion(), libMesh::AdvectionSystem::assemble_claw_rhs(), libMesh::FEMSystem::assemble_qoi(), libMesh::Nemesis_IO::assert_symmetric_cmaps(), libMesh::MeshCommunication::assign_global_indices(), libMesh::Partitioner::assign_partitioning(), libMesh::MeshTools::Generation::build_extrusion(), libMesh::Partitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::PetscDMWrapper::build_section(), libMesh::PetscDMWrapper::build_sf(), libMesh::MeshBase::cache_elem_data(), libMesh::DofMap::check_dirichlet_bcid_consistency(), libMesh::MeshTetInterface::check_hull_integrity(), libMesh::MeshBase::complete_preparation(), libMesh::RBConstruction::compute_Fq_representor_innerprods(), libMesh::RBConstruction::compute_max_error_bound(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::RBConstruction::compute_output_dual_innerprods(), libMesh::RBConstruction::compute_residual_dual_norm_slow(), libMesh::RBSCMConstruction::compute_SCM_bounds_on_training_set(), libMesh::DofMap::computed_sparsity_already(), libMesh::Problem_Interface::computeJacobian(), libMesh::Problem_Interface::computePreconditioner(), PetscSolverConfiguration::configure_solver(), libMesh::ContinuationSystem::ContinuationSystem(), libMesh::MeshBase::copy_constraint_rows(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::ExodusII_IO::copy_nodal_solution(), libMesh::ExodusII_IO::copy_scalar_solution(), libMesh::CondensedEigenSystem::copy_super_to_sub(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::MeshTools::create_bounding_box(), libMesh::DofMap::create_dof_constraints(), libMesh::MeshTools::create_nodal_bounding_box(), libMesh::MeshRefinement::create_parent_error_vector(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::MeshTools::create_subdomain_bounding_box(), libMesh::PetscMatrix< T >::create_submatrix_nosort(), create_wrapped_function(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::MeshBase::detect_interior_parents(), libMesh::RBEIMEvaluation::distribute_bfs(), DMlibMeshFunction(), DMlibMeshJacobian(), DMlibMeshSetSystem_libMesh(), DMVariableBounds_libMesh(), libMesh::DTKSolutionTransfer::DTKSolutionTransfer(), libMesh::MeshRefinement::eliminate_unrefined_patches(), libMesh::RBEIMConstruction::enrich_eim_approximation_on_interiors(), libMesh::RBEIMConstruction::enrich_eim_approximation_on_nodes(), libMesh::RBEIMConstruction::enrich_eim_approximation_on_sides(), libMesh::TransientRBConstruction::enrich_RB_space(), libMesh::EpetraVector< T >::EpetraVector(), AssembleOptimization::equality_constraints(), libMesh::AdjointRefinementEstimator::estimate_error(), libMesh::ExactErrorEstimator::estimate_error(), libMesh::JumpErrorEstimator::estimate_error(), libMesh::PatchRecoveryErrorEstimator::estimate_error(), libMesh::WeightedPatchRecoveryErrorEstimator::estimate_error(), libMesh::SmoothnessEstimator::estimate_smoothness(), libMesh::MeshRefinement::flag_elements_by_elem_fraction(), libMesh::MeshRefinement::flag_elements_by_error_fraction(), libMesh::MeshRefinement::flag_elements_by_error_tolerance(), libMesh::MeshRefinement::flag_elements_by_mean_stddev(), libMesh::MeshRefinement::flag_elements_by_nelem_target(), libMesh::RBEIMEvaluation::gather_bfs(), libMesh::DofMap::gather_constraints(), libMesh::MeshfreeInterpolation::gather_remote_data(), libMesh::CondensedEigenSystem::get_eigenpair(), libMesh::RBEIMEvaluation::get_eim_basis_function_node_value(), libMesh::RBEIMEvaluation::get_eim_basis_function_side_value(), libMesh::RBEIMEvaluation::get_eim_basis_function_value(), libMesh::MeshBase::get_info(), libMesh::RBEIMEvaluation::get_interior_basis_functions_as_vecs(), libMesh::ImplicitSystem::get_linear_solver(), libMesh::RBEIMConstruction::get_max_abs_value(), libMesh::RBEIMConstruction::get_node_max_abs_value(), libMesh::RBEIMEvaluation::get_parametrized_function_node_value(), libMesh::RBEIMEvaluation::get_parametrized_function_side_value(), libMesh::RBEIMEvaluation::get_parametrized_function_value(), libMesh::RBEIMConstruction::get_random_point(), libMesh::RBEIMConstruction::get_random_point(), libMesh::RBEIMConstruction::get_random_point(), libMesh::MeshTetInterface::improve_hull_integrity(), AssembleOptimization::inequality_constraints(), AssembleOptimization::inequality_constraints_jacobian(), libMesh::StaticCondensation::init(), libMesh::TimeSolver::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::LocationMap< T >::init(), libMesh::PetscDMWrapper::init_and_attach_petscdm(), libMesh::PetscDMWrapper::init_and_attach_petscdm(), ElasticitySystem::init_data(), libMesh::AdvectionSystem::init_data(), libMesh::ClawSystem::init_data(), libMesh::PetscDMWrapper::init_petscdm(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::OptimizationSystem::initialize_equality_constraints_storage(), libMesh::OptimizationSystem::initialize_inequality_constraints_storage(), libMesh::RBEIMConstruction::initialize_parametrized_functions_in_training_set(), libMesh::RBEIMConstruction::inner_product(), integrate_function(), libMesh::MeshTools::Modification::interpolate_surface(), libMesh::MeshTools::libmesh_assert_consistent_distributed(), libMesh::MeshTools::libmesh_assert_consistent_distributed_nodes(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_equal_connectivity(), libMesh::MeshTools::libmesh_assert_equal_points(), libMesh::MeshTools::libmesh_assert_parallel_consistent_new_node_procids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_constraint_rows(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_flags(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_p_levels(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshTools::libmesh_assert_valid_unique_ids(), libMesh::libmesh_petsc_linesearch_shellfunc(), libMesh::libmesh_petsc_preconditioner_apply(), libMesh::libmesh_petsc_snes_mffd_interface(), libMesh::libmesh_petsc_snes_postcheck(), libMesh::MeshRefinement::limit_level_mismatch_at_edge(), libMesh::MeshRefinement::limit_level_mismatch_at_node(), libMesh::MeshRefinement::limit_overrefined_boundary(), libMesh::MeshRefinement::limit_underrefined_boundary(), libMesh::LinearImplicitSystem::LinearImplicitSystem(), main(), libMesh::MeshCommunication::make_elems_parallel_consistent(), libMesh::MeshCommunication::make_new_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_new_nodes_parallel_consistent(), libMesh::MeshCommunication::make_node_bcids_parallel_consistent(), libMesh::MeshCommunication::make_node_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_proc_ids_parallel_consistent(), libMesh::MeshCommunication::make_node_unique_ids_parallel_consistent(), libMesh::MeshCommunication::make_nodes_parallel_consistent(), libMesh::MeshCommunication::make_p_levels_parallel_consistent(), libMesh::TransientRBConstruction::mass_matrix_scaled_matvec(), libMesh::FEMSystem::mesh_position_set(), libMesh::TriangulatorInterface::MeshedHole::MeshedHole(), LinearElasticityWithContact::move_mesh(), libMesh::DistributedMesh::n_active_elem(), libMesh::MeshTools::n_active_levels(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshTools::n_connected_components(), libMesh::DofMap::n_constrained_dofs(), libMesh::MeshBase::n_constraint_rows(), libMesh::DofMap::n_dofs(), libMesh::DofMap::n_dofs_per_processor(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::CondensedEigenSystem::n_global_non_condensed_dofs(), libMesh::MeshTools::n_levels(), MixedOrderTest::n_neighbor_links(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::SparsityPattern::Build::n_nonzeros(), libMesh::MeshTools::n_p_levels(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::RBEIMEvaluation::node_distribute_bfs(), libMesh::RBEIMEvaluation::node_gather_bfs(), libMesh::RBEIMConstruction::node_inner_product(), libMesh::PetscVector< T >::operator=(), libMesh::MeshBase::operator==(), libMesh::DistributedMesh::parallel_max_elem_id(), libMesh::DistributedMesh::parallel_max_node_id(), libMesh::DistributedMesh::parallel_max_unique_id(), libMesh::ReplicatedMesh::parallel_max_unique_id(), libMesh::DistributedMesh::parallel_n_elem(), libMesh::DistributedMesh::parallel_n_nodes(), libMesh::SparsityPattern::Build::parallel_sync(), libMesh::BoundaryInfo::parallel_sync_node_ids(), libMesh::BoundaryInfo::parallel_sync_side_ids(), libMesh::MeshTools::paranoid_n_levels(), libMesh::Partitioner::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::petsc_auto_fieldsplit(), LaplaceSystem::postprocess(), PoissonSystem::postprocess(), libMesh::MeshBase::print_constraint_rows(), libMesh::DofMap::print_dof_constraints(), libMesh::DofMap::process_mesh_constraint_rows(), libMesh::Partitioner::processor_pairs_to_interface_nodes(), libMesh::InterMeshProjection::project_system_vectors(), libMesh::XdrIO::read(), libMesh::Nemesis_IO::read(), FEMParameters::read(), libMesh::EquationSystems::read(), libMesh::CheckpointIO::read_header(), libMesh::ExodusII_IO::read_header(), libMesh::System::read_header(), libMesh::XdrIO::read_header(), libMesh::RBEIMEvaluation::read_in_interior_basis_functions(), libMesh::RBEIMEvaluation::read_in_node_basis_functions(), libMesh::RBEIMEvaluation::read_in_side_basis_functions(), libMesh::RBEvaluation::read_in_vectors_from_multiple_files(), libMesh::RBConstruction::read_riesz_representors_from_files(), libMesh::TransientRBConstruction::read_riesz_representors_from_files(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::Nemesis_IO_Helper::read_var_names_impl(), MeshFunctionTest::read_variable_info_from_output_data(), libMesh::MeshBase::recalculate_n_partitions(), libMesh::SimplexRefiner::refine_via_edges(), libMesh::StaticCondensationDofMap::reinit(), libMesh::BoundaryInfo::remove_edge_id(), libMesh::BoundaryInfo::remove_node_id(), libMesh::BoundaryInfo::remove_shellface_id(), libMesh::BoundaryInfo::remove_side_id(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::DistributedMesh::renumber_nodes_and_elements(), LinearElasticityWithContact::residual_and_jacobian(), OverlappingAlgebraicGhostingTest::run_ghosting_test(), OverlappingCouplingGhostingTest::run_sparsity_pattern_test(), scale_mesh_and_plot(), libMesh::DofMap::scatter_constraints(), libMesh::CheckpointIO::select_split_config(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::send_and_insert_dof_values(), libMesh::TransientRBConstruction::set_error_temporal_data(), libMesh::Partitioner::set_interface_node_processor_ids_BFS(), libMesh::Partitioner::set_interface_node_processor_ids_linear(), libMesh::Partitioner::set_interface_node_processor_ids_petscpartitioner(), libMesh::Partitioner::set_node_processor_ids(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::Partitioner::set_parent_processor_ids(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::PetscDiffSolver::setup_petsc_data(), libMesh::RBEIMEvaluation::side_distribute_bfs(), libMesh::RBEIMEvaluation::side_gather_bfs(), libMesh::RBEIMConstruction::side_inner_product(), libMesh::Partitioner::single_partition(), libMesh::LaplaceMeshSmoother::smooth(), libMesh::VariationalMeshSmoother::smooth(), libMesh::NoxNonlinearSolver< Number >::solve(), libMesh::ClawSystem::solve_conservation_law(), libMesh::split_mesh(), libMesh::RBEIMConstruction::store_eim_solutions_for_training_set(), libMesh::MeshBase::subdomain_ids(), libMesh::BoundaryInfo::sync(), libMesh::MeshBase::sync_subdomain_name_map(), ConstraintOperatorTest::test1DCoarseningNewNodes(), ConstraintOperatorTest::test1DCoarseningOperator(), MeshFunctionTest::test_bad_gradient_var_with_out_of_mesh_value(), MeshFunctionTest::test_bad_hessian_var_with_out_of_mesh_value(), MeshfunctionDFEM::test_mesh_function_dfem(), MeshfunctionDFEM::test_mesh_function_dfem_grad(), MeshFunctionTest::test_p_level(), ExodusC0PolyhedronTest::test_write_and_read_hexagonal_prism(), ExodusC0PolygonTest::test_write_and_read_pentagon(), DofMapTest::testBadElemFECombo(), SystemsTest::testBlockRestrictedVarNDofs(), BoundaryInfoTest::testBoundaryOnChildrenErrors(), CheckpointIOTest::testC0PolygonCheckpoint(), VolumeTest::testC0PolygonMethods(), CheckpointIOTest::testC0PolyhedronCheckpoint(), VolumeTest::testC0PolyhedronMethods(), ConstraintOperatorTest::testCoreform(), ConnectedComponentsTest::testEdge(), MeshInputTest::testExodusIGASidesets(), MeshTriangulationTest::testFoundCenters(), PointLocatorTest::testLocator(), BoundaryInfoTest::testMesh(), BoundaryMeshSubdomainTest::testPerBoundarySubdomain(), PointLocatorTest::testPlanar(), MeshTriangulationTest::testPoly2TriEdge3ToTri7CenterFixup(), MeshTriangulationTest::testPoly2TriRefinementBase(), SystemsTest::testProjectCubeWithMeshFunction(), SystemsTest::testProjectScalarCoarsening(), BoundaryInfoTest::testRenumber(), BoundaryInfoTest::testSelectiveRenumber(), BoundaryMeshSubdomainTest::testSingleSubdomain(), CheckpointIOTest::testSplitter(), MeshInputTest::testTetgenIO(), MeshTriangulationTest::testTriangulatorInterp(), MeshTriangulationTest::testTriangulatorMeshedHoles(), MeshTriangulationTest::testTriangulatorRoundHole(), MeshSmootherTest::testVariationalSmoother(), libMesh::MeshTools::total_weight(), libMesh::RBConstruction::train_reduced_basis_with_POD(), libMesh::MeshfreeSolutionTransfer::transfer(), libMesh::MeshFunctionSolutionTransfer::transfer(), libMesh::Poly2TriTriangulator::triangulate(), libMesh::RBConstruction::truth_assembly(), libMesh::TransientRBConstruction::truth_assembly(), update_current_local_solution(), libMesh::TransientRBConstruction::update_RB_initial_condition_all_N(), libMesh::RBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_residual_terms(), libMesh::RBConstruction::update_residual_terms(), libMesh::MeshTools::volume(), libMesh::STLIO::write(), libMesh::XdrIO::write(), libMesh::NameBasedIO::write(), libMesh::VTKIO::write_nodal_data(), libMesh::RBEIMEvaluation::write_out_interior_basis_functions(), libMesh::RBEIMEvaluation::write_out_node_basis_functions(), libMesh::RBEIMEvaluation::write_out_side_basis_functions(), libMesh::RBEvaluation::write_out_vectors(), libMesh::RBConstruction::write_riesz_representors_to_files(), libMesh::TransientRBConstruction::write_riesz_representors_to_files(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::RBDataSerialization::RBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::TransientRBEvaluationSerialization::write_to_file(), libMesh::RBDataSerialization::RBEIMEvaluationSerialization::write_to_file(), and libMesh::RBDataSerialization::RBSCMEvaluationSerialization::write_to_file().

◆ copy_from_hash()

template<typename T >
std::unique_ptr< PetscMatrix< T > > libMesh::PetscMatrix< T >::copy_from_hash ( )

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

Definition at line 1254 of file petsc_matrix.C.

1255{
1256 Mat xaij;
1257 libmesh_assert(this->initialized());
1258 libmesh_assert(!this->closed());
1259 LibmeshPetscCall(MatDuplicate(this->_mat, MAT_DO_NOT_COPY_VALUES, &xaij));
1260 LibmeshPetscCall(MatCopyHashToXAIJ(this->_mat, xaij));
1261 return std::make_unique<PetscMatrix<T>>(xaij, this->comm(), /*destroy_on_exit=*/true);
1262}

References libMesh::closed(), libMesh::initialized(), and libMesh::libmesh_assert().

Referenced by PetscMatrixTest::testPetscCopyFromHash().

◆ create_submatrix()

template<typename T >
virtual void libMesh::SparseMatrix< T >::create_submatrix ( SparseMatrix< T > &  submatrix,
const std::vector< numeric_index_type > &  rows,
const std::vector< numeric_index_type > &  cols 
) const
inlinevirtualinherited

This function creates a matrix called "submatrix" which is defined by the row and column indices given in the "rows" and "cols" entries.

Currently this operation is only defined for the PetscMatrixBase subclasses. Note: The rows and cols vectors need to be sorted; Use the nosort version below if rows and cols vectors are not sorted; The rows and cols only contain indices that are owned by this processor.

Definition at line 534 of file sparse_matrix.h.

537 {
538 this->_get_submatrix(submatrix,
539 rows,
540 cols,
541 false); // false means DO NOT REUSE submatrix
542 }
virtual void _get_submatrix(SparseMatrix< T > &, const std::vector< numeric_index_type > &, const std::vector< numeric_index_type > &, const bool) const
Protected implementation of the create_submatrix and reinit_submatrix routines.

References libMesh::SparseMatrix< T >::_get_submatrix().

Referenced by libMesh::CondensedEigenSystem::copy_super_to_sub(), libMesh::libmesh_petsc_DMCreateInterpolation(), and libMesh::CondensedEigenSystem::solve().

◆ create_submatrix_nosort()

template<typename T >
void libMesh::PetscMatrix< T >::create_submatrix_nosort ( SparseMatrix< T > &  submatrix,
const std::vector< numeric_index_type > &  rows,
const std::vector< numeric_index_type > &  cols 
) const
overridevirtual

Similar to the create_submatrix function, this function creates a submatrix which is defined by the indices given in the rows and cols vectors.

Note: Both rows and cols can be unsorted; Use the above function for better efficiency if your indices are sorted; rows and cols can contain indices that are owned by other processors.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 853 of file petsc_matrix.C.

856{
857 if (!this->closed())
858 {
859 libmesh_deprecated();
860 libmesh_warning("The matrix must be assembled before calling PetscMatrix::create_submatrix_nosort().\n"
861 "Please update your code, as this warning will become an error in a future release.");
862 const_cast<PetscMatrix<T> *>(this)->close();
863 }
864
865 // Make sure the SparseMatrix passed in is really a PetscMatrix
866 PetscMatrix<T> * petsc_submatrix = cast_ptr<PetscMatrix<T> *>(&submatrix);
867
868 LibmeshPetscCall(MatZeroEntries(petsc_submatrix->_mat));
869
870 PetscInt pc_ncols = 0;
871 const PetscInt * pc_cols;
872 const PetscScalar * pc_vals;
873
874 // // data for creating the submatrix
875 std::vector<PetscInt> sub_cols;
876 std::vector<PetscScalar> sub_vals;
877
878 for (auto i : index_range(rows))
879 {
880 PetscInt sub_rid[] = {static_cast<PetscInt>(i)};
881 PetscInt rid = static_cast<PetscInt>(rows[i]);
882 // only get value from local rows, and set values to the corresponding columns in the submatrix
883 if (rows[i]>= this->row_start() && rows[i]< this->row_stop())
884 {
885 // get one row of data from the original matrix
886 LibmeshPetscCall(MatGetRow(this->_mat, rid, &pc_ncols, &pc_cols, &pc_vals));
887 // extract data from certain cols, save the indices and entries sub_cols and sub_vals
888 for (auto j : index_range(cols))
889 {
890 for (unsigned int idx = 0; idx< static_cast<unsigned int>(pc_ncols); idx++)
891 {
892 if (pc_cols[idx] == static_cast<PetscInt>(cols[j]))
893 {
894 sub_cols.push_back(static_cast<PetscInt>(j));
895 sub_vals.push_back(pc_vals[idx]);
896 }
897 }
898 }
899 // set values
900 LibmeshPetscCall(MatSetValues(petsc_submatrix->_mat,
901 1,
902 sub_rid,
903 static_cast<PetscInt>(sub_vals.size()),
904 sub_cols.data(),
905 sub_vals.data(),
906 INSERT_VALUES));
907 LibmeshPetscCall(MatRestoreRow(this->_mat, rid, &pc_ncols, &pc_cols, &pc_vals));
908 // clear data for this row
909 sub_cols.clear();
910 sub_vals.clear();
911 }
912 }
913 MatAssemblyBeginEnd(petsc_submatrix->comm(), petsc_submatrix->_mat, MAT_FINAL_ASSEMBLY);
914 // Specify that the new submatrix is initialized and close it.
915 petsc_submatrix->_is_initialized = true;
916 petsc_submatrix->close();
917}
virtual numeric_index_type row_start() const override
virtual numeric_index_type row_stop() const override
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.

References libMesh::SparseMatrix< T >::_is_initialized, libMesh::PetscMatrixBase< T >::_mat, libMesh::PetscMatrixBase< T >::close(), libMesh::closed(), libMesh::ParallelObject::comm(), and libMesh::index_range().

◆ disable_print_counter_info()

void libMesh::ReferenceCounter::disable_print_counter_info ( )
staticinherited

Definition at line 100 of file reference_counter.C.

101{
102 _enable_print_counter = false;
103 return;
104}
static bool _enable_print_counter
Flag to control whether reference count information is printed when print_info is called.

References libMesh::ReferenceCounter::_enable_print_counter.

◆ enable_print_counter_info()

void libMesh::ReferenceCounter::enable_print_counter_info ( )
staticinherited

Methods to enable/disable the reference counter output from print_info().

Enabled by default.

Definition at line 94 of file reference_counter.C.

95{
97 return;
98}

References libMesh::ReferenceCounter::_enable_print_counter.

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ finish_initialization()

template<typename T >
void libMesh::PetscMatrix< T >::finish_initialization ( )

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

Definition at line 329 of file petsc_matrix.C.

330{
331 // Make it an error for PETSc to allocate new nonzero entries during assembly. For old PETSc
332 // versions this option must be set after preallocation for MPIAIJ matrices
333 LibmeshPetscCall(MatSetOption(this->_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE));
334 this->_is_initialized = true;
335}

Referenced by PetscMatrixTest::testPetscCopyFromHash().

◆ flush()

template<typename T >
void libMesh::PetscMatrix< T >::flush ( )
overridevirtual

For PETSc matrix , this function is similar to close but without shrinking memory.

This is useful when we want to switch between ADD_VALUES and INSERT_VALUES. close should be called before using the matrix.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 962 of file petsc_matrix.C.

963{
964 semiparallel_only();
965
966 MatAssemblyBeginEnd(this->comm(), this->_mat, MAT_FLUSH_ASSEMBLY);
967}

◆ frobenius_norm()

template<typename T >
Real libMesh::PetscMatrix< T >::frobenius_norm ( ) const

Definition at line 497 of file petsc_matrix.C.

498{
500}

◆ get_context()

template<typename T >
PetscMatrixBase< T > * libMesh::PetscMatrixBase< T >::get_context ( Mat  mat,
const TIMPI::Communicator comm 
)
staticinherited
Returns
The context for mat if it exists, else a nullptr

Definition at line 116 of file petsc_matrix_base.C.

117{
118 void * ctx;
119 PetscContainer container;
120 LibmeshPetscCall2(comm, PetscObjectQuery((PetscObject)mat, "PetscMatrixCtx", (PetscObject *)&container));
121 if (!container)
122 return nullptr;
123
124 LibmeshPetscCall2(comm, PetscContainerGetPointer(container, &ctx));
126 return static_cast<PetscMatrixBase<T> *>(ctx);
127}

References libMesh::ctx, and libMesh::libmesh_assert().

Referenced by DMlibMeshJacobian(), and form_matrixA().

◆ get_diagonal()

template<typename T >
void libMesh::PetscMatrix< T >::get_diagonal ( NumericVector< T > &  dest) const
overridevirtual

Copies the diagonal part of the matrix into dest.

Implements libMesh::SparseMatrix< T >.

Definition at line 921 of file petsc_matrix.C.

922{
923 // Make sure the NumericVector passed in is really a PetscVector
924 PetscVector<T> & petsc_dest = cast_ref<PetscVector<T> &>(dest);
925
926 // Needs a const_cast since PETSc does not work with const.
927 LibmeshPetscCall(MatGetDiagonal(const_cast<PetscMatrix<T> *>(this)->mat(),petsc_dest.vec()));
928}
This class provides a nice interface to PETSc's Vec object.

References libMesh::PetscVector< T >::vec().

◆ get_info()

std::string libMesh::ReferenceCounter::get_info ( )
staticinherited

Gets a string containing the reference information.

Definition at line 47 of file reference_counter.C.

48{
49#if defined(LIBMESH_ENABLE_REFERENCE_COUNTING) && defined(DEBUG)
50
51 std::ostringstream oss;
52
53 oss << '\n'
54 << " ---------------------------------------------------------------------------- \n"
55 << "| Reference count information |\n"
56 << " ---------------------------------------------------------------------------- \n";
57
58 for (const auto & [name, cd] : _counts)
59 oss << "| " << name << " reference count information:\n"
60 << "| Creations: " << cd.first << '\n'
61 << "| Destructions: " << cd.second << '\n';
62
63 oss << " ---------------------------------------------------------------------------- \n";
64
65 return oss.str();
66
67#else
68
69 return "";
70
71#endif
72}
static Counts _counts
Actually holds the data.
std::string name(const ElemQuality q)
This function returns a string containing some name for q.

References libMesh::ReferenceCounter::_counts.

Referenced by libMesh::ReferenceCounter::print_info().

◆ get_row()

template<typename T >
void libMesh::PetscMatrix< T >::get_row ( numeric_index_type  i,
std::vector< numeric_index_type > &  indices,
std::vector< T > &  values 
) const
overridevirtual

Get a row from the matrix.

Parameters
iThe matrix row to get
indicesA container that will be filled with the column indices corresponding to (possibly) non-zero values
valuesA container holding the column values

Implements libMesh::SparseMatrix< T >.

Definition at line 1148 of file petsc_matrix.C.

1151{
1152 libmesh_assert (this->initialized());
1153
1154 const PetscScalar * petsc_row;
1155 const PetscInt * petsc_cols;
1156
1157 PetscInt
1158 ncols=0,
1159 i_val = static_cast<PetscInt>(i_in);
1160
1161 // the matrix needs to be closed for this to work
1162 // this->close();
1163 // but closing it is a semiparallel operation; we want operator()
1164 // to run on one processor.
1165 libmesh_assert(this->closed());
1166
1167 // PETSc makes no effort at being thread safe. Helgrind complains about
1168 // possible data races even just in PetscFunctionBegin (due to things
1169 // like stack counter incrementing). Perhaps we could ignore
1170 // this, but there are legitimate data races for Mat data members like
1171 // mat->getrowactive between MatGetRow and MatRestoreRow. Moreover,
1172 // there could be a write into mat->rowvalues during MatGetRow from
1173 // one thread while we are attempting to read from mat->rowvalues
1174 // (through petsc_cols) during data copy in another thread. So
1175 // the safe thing to do is to lock the whole method
1176
1177 std::lock_guard<std::mutex> lock(_petsc_matrix_mutex);
1178
1179 LibmeshPetscCall(MatGetRow(this->_mat, i_val, &ncols, &petsc_cols, &petsc_row));
1180
1181 // Copy the data
1182 indices.resize(static_cast<std::size_t>(ncols));
1183 values.resize(static_cast<std::size_t>(ncols));
1184
1185 for (auto i : index_range(indices))
1186 {
1187 indices[i] = static_cast<numeric_index_type>(petsc_cols[i]);
1188 values[i] = static_cast<T>(petsc_row[i]);
1189 }
1190
1191 LibmeshPetscCall(MatRestoreRow(this->_mat, i_val,
1192 &ncols, &petsc_cols, &petsc_row));
1193}

References libMesh::closed(), libMesh::index_range(), libMesh::initialized(), and libMesh::libmesh_assert().

◆ get_transpose()

template<typename T >
void libMesh::PetscMatrix< T >::get_transpose ( SparseMatrix< T > &  dest) const
overridevirtual

Copies the transpose of the matrix into dest, which may be *this.

Implements libMesh::SparseMatrix< T >.

Definition at line 933 of file petsc_matrix.C.

934{
935 // Make sure the SparseMatrix passed in is really a PetscMatrix
936 PetscMatrix<T> & petsc_dest = cast_ref<PetscMatrix<T> &>(dest);
937
938 // If we aren't reusing the matrix then need to clear dest,
939 // otherwise we get a memory leak
940 if (&petsc_dest != this)
941 dest.clear();
942
943 if (&petsc_dest == this)
944 // The MAT_REUSE_MATRIX flag was replaced by MAT_INPLACE_MATRIX
945 // in PETSc 3.7.0
946#if PETSC_VERSION_LESS_THAN(3,7,0)
947 LibmeshPetscCall(MatTranspose(this->_mat,MAT_REUSE_MATRIX, &petsc_dest._mat));
948#else
949 LibmeshPetscCall(MatTranspose(this->_mat, MAT_INPLACE_MATRIX, &petsc_dest._mat));
950#endif
951 else
952 LibmeshPetscCall(MatTranspose(this->_mat,MAT_INITIAL_MATRIX, &petsc_dest._mat));
953
954 // Specify that the transposed matrix is initialized and close it.
955 petsc_dest._is_initialized = true;
956 petsc_dest.close();
957}

References libMesh::SparseMatrix< T >::_is_initialized, libMesh::PetscMatrixBase< T >::_mat, libMesh::SparseMatrix< T >::clear(), and libMesh::PetscMatrixBase< T >::close().

◆ increment_constructor_count()

void libMesh::ReferenceCounter::increment_constructor_count ( const std::string &  name)
inlineprotectednoexceptinherited

Increments the construction counter.

Should be called in the constructor of any derived class that will be reference counted.

Definition at line 183 of file reference_counter.h.

184{
185 libmesh_try
186 {
187 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
188 std::pair<unsigned int, unsigned int> & p = _counts[name];
189 p.first++;
190 }
191 libmesh_catch (...)
192 {
193 auto stream = libMesh::err.get();
194 stream->exceptions(stream->goodbit); // stream must not throw
195 libMesh::err << "Encountered unrecoverable error while calling "
196 << "ReferenceCounter::increment_constructor_count() "
197 << "for a(n) " << name << " object." << std::endl;
198 std::terminate();
199 }
200}
streamT * get()
Rather than implement every ostream/ios/ios_base function, we'll be lazy and make esoteric uses go th...
spin_mutex spin_mtx
A convenient spin mutex object which can be used for obtaining locks.
Definition threads.C:30
OStreamProxy err

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), libMesh::ReferenceCountedObject< T >::ReferenceCountedObject(), and libMesh::ReferenceCountedObject< T >::ReferenceCountedObject().

◆ increment_destructor_count()

void libMesh::ReferenceCounter::increment_destructor_count ( const std::string &  name)
inlineprotectednoexceptinherited

Increments the destruction counter.

Should be called in the destructor of any derived class that will be reference counted.

Definition at line 207 of file reference_counter.h.

208{
209 libmesh_try
210 {
211 Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
212 std::pair<unsigned int, unsigned int> & p = _counts[name];
213 p.second++;
214 }
215 libmesh_catch (...)
216 {
217 auto stream = libMesh::err.get();
218 stream->exceptions(stream->goodbit); // stream must not throw
219 libMesh::err << "Encountered unrecoverable error while calling "
220 << "ReferenceCounter::increment_destructor_count() "
221 << "for a(n) " << name << " object." << std::endl;
222 std::terminate();
223 }
224}

References libMesh::err, libMesh::BasicOStreamProxy< charT, traits >::get(), and libMesh::Threads::spin_mtx.

Referenced by libMesh::ReferenceCountedObject< T >::~ReferenceCountedObject().

◆ init() [1/3]

template<typename T >
void libMesh::PetscMatrix< T >::init ( const numeric_index_type  m,
const numeric_index_type  n,
const numeric_index_type  m_l,
const numeric_index_type  n_l,
const numeric_index_type  n_nz = 30,
const numeric_index_type  n_oz = 10,
const numeric_index_type  blocksize = 1 
)
overridevirtual

Initialize a PETSc matrix.

Parameters
mThe global number of rows.
nThe global number of columns.
m_lThe local number of rows.
n_lThe local number of columns.
n_nzThe number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
n_ozThe number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
blocksizeOptional value indicating dense coupled blocks for systems with multiple variables all of the same type.

Implements libMesh::SparseMatrix< T >.

Definition at line 214 of file petsc_matrix.C.

221{
222 this->init_without_preallocation(m_in, n_in, m_l, n_l, blocksize_in);
223
224 PetscInt n_nz = static_cast<PetscInt>(nnz);
225 PetscInt n_oz = static_cast<PetscInt>(noz);
226
227#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
228 if (blocksize > 1)
229 {
230 LibmeshPetscCall(MatSeqBAIJSetPreallocation(this->_mat, blocksize, n_nz/blocksize, NULL));
231 LibmeshPetscCall(MatMPIBAIJSetPreallocation(this->_mat, blocksize,
232 n_nz/blocksize, NULL,
233 n_oz/blocksize, NULL));
234 }
235 else
236#endif
237 {
238 switch (this->_mat_type) {
239 case AIJ:
240 LibmeshPetscCall(MatSeqAIJSetPreallocation(this->_mat, n_nz, NULL));
241 LibmeshPetscCall(MatMPIAIJSetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
242 break;
243
244 case HYPRE:
245#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
246 LibmeshPetscCall(MatHYPRESetPreallocation(this->_mat, n_nz, NULL, n_oz, NULL));
247#else
248 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
249#endif
250 break;
251
252 default: libmesh_error_msg("Unsupported petsc matrix type");
253 }
254 }
255
256 this->finish_initialization();
257}
void finish_initialization()
Finish up the initialization process.
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)
Perform matrix initialization steps sans preallocation.

References libMesh::AIJ, finish_initialization(), and libMesh::HYPRE.

Referenced by libMesh::PetscMatrix< T >::PetscMatrix().

◆ init() [2/3]

template<typename T >
void libMesh::PetscMatrix< T >::init ( const numeric_index_type  m,
const numeric_index_type  n,
const numeric_index_type  m_l,
const numeric_index_type  n_l,
const std::vector< numeric_index_type > &  n_nz,
const std::vector< numeric_index_type > &  n_oz,
const numeric_index_type  blocksize = 1 
)

Initialize a PETSc matrix.

Parameters
mThe global number of rows.
nThe global number of columns.
m_lThe local number of rows.
n_lThe local number of columns.
n_nzarray containing the number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
n_ozArray containing the number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
blocksizeOptional value indicating dense coupled blocks for systems with multiple variables all of the same type.

Definition at line 338 of file petsc_matrix.C.

345{
346 this->init_without_preallocation(m_in, n_in, m_l, n_l, blocksize_in);
347 this->preallocate(m_l, n_nz, n_oz, blocksize_in);
348
349 this->finish_initialization();
350}
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)

References finish_initialization().

◆ init() [3/3]

template<typename T >
virtual void libMesh::PetscMatrix< T >::init ( ParallelType  type = PARALLEL)
overridevirtual

Initialize this matrix using the sparsity structure computed by dof_map.

Parameters
typeThe serial/parallel/ghosted type of the matrix

Implements libMesh::SparseMatrix< T >.

◆ init_without_preallocation()

template<typename T >
void libMesh::PetscMatrix< T >::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 
)

Perform matrix initialization steps sans preallocation.

Parameters
mThe global number of rows
nThe global number of columns
m_lThe local number of rows
n_lThe local number of columns
blocksizeThe matrix block size

Definition at line 135 of file petsc_matrix.C.

140{
141 // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
142 libmesh_ignore(blocksize_in);
143
144 // Clear initialized matrices
145 if (this->initialized())
146 this->clear();
147
148 PetscInt m_global = static_cast<PetscInt>(m_in);
149 PetscInt n_global = static_cast<PetscInt>(n_in);
150 PetscInt m_local = static_cast<PetscInt>(m_l);
151 PetscInt n_local = static_cast<PetscInt>(n_l);
152
153 LibmeshPetscCall(MatCreate(this->comm().get(), &this->_mat));
154 LibmeshPetscCall(MatSetSizes(this->_mat, m_local, n_local, m_global, n_global));
155 PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
156 LibmeshPetscCall(MatSetBlockSize(this->_mat,blocksize));
157 LibmeshPetscCall(MatSetOptionsPrefix(this->_mat, ""));
158
159#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
160 if (blocksize > 1)
161 {
162 // specified blocksize, bs>1.
163 // double check sizes.
164 libmesh_assert_equal_to (m_local % blocksize, 0);
165 libmesh_assert_equal_to (n_local % blocksize, 0);
166 libmesh_assert_equal_to (m_global % blocksize, 0);
167 libmesh_assert_equal_to (n_global % blocksize, 0);
168 libmesh_assert_equal_to (n_nz % blocksize, 0);
169 libmesh_assert_equal_to (n_oz % blocksize, 0);
170
171 LibmeshPetscCall(MatSetType(this->_mat, MATBAIJ)); // Automatically chooses seqbaij or mpibaij
172
173 // MatSetFromOptions needs to happen before Preallocation routines
174 // since MatSetFromOptions can change matrix type and remove incompatible
175 // preallocation
176 LibmeshPetscCall(MatSetFromOptions(this->_mat));
177 }
178 else
179#endif
180 {
181 switch (this->_mat_type) {
182 case AIJ:
183 LibmeshPetscCall(MatSetType(this->_mat, MATAIJ)); // Automatically chooses seqaij or mpiaij
184
185 // MatSetFromOptions needs to happen before Preallocation routines
186 // since MatSetFromOptions can change matrix type and remove incompatible
187 // preallocation
188 LibmeshPetscCall(MatSetFromOptions(this->_mat));
189 break;
190
191 case HYPRE:
192#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
193 LibmeshPetscCall(MatSetType(this->_mat, MATHYPRE));
194
195 // MatSetFromOptions needs to happen before Preallocation routines
196 // since MatSetFromOptions can change matrix type and remove incompatible
197 // preallocation
198 LibmeshPetscCall(MatSetFromOptions(this->_mat));
199#else
200 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
201#endif
202 break;
203
204 default: libmesh_error_msg("Unsupported petsc matrix type");
205 }
206 }
207
208 this->set_context ();
209}
void set_context()
Set the context (ourself) for _mat.
virtual void clear() noexcept override
clear() is called from the destructor, so it should not throw.

References libMesh::AIJ, libMesh::HYPRE, libMesh::initialized(), and libMesh::libmesh_ignore().

Referenced by PetscMatrixTest::testPetscCopyFromHash().

◆ initialized()

template<typename T >
virtual bool libMesh::SparseMatrix< T >::initialized ( ) const
inlinevirtualinherited
Returns
true if the matrix has been initialized, false otherwise.

Reimplemented in libMesh::StaticCondensation.

Definition at line 133 of file sparse_matrix.h.

133{ return _is_initialized; }

References libMesh::SparseMatrix< T >::_is_initialized.

Referenced by libMesh::PetscMatrix< T >::_get_submatrix(), libMesh::ImplicitSystem::assemble(), libMesh::System::init_matrices(), and libMesh::StaticCondensation::initialized().

◆ l1_norm()

template<typename T >
Real libMesh::PetscMatrix< T >::l1_norm ( ) const
overridevirtual
Returns
The \( \ell_1 \)-norm of the matrix, that is the max column sum: \( |M|_1 = \max_{j} \sum_{i} |M_{ij}| \)

This is the natural matrix norm that is compatible with the \( \ell_1 \)-norm for vectors, i.e. \( |Mv|_1 \leq |M|_1 |v|_1 \). (cf. Haemmerlin-Hoffmann : Numerische Mathematik)

Implements libMesh::SparseMatrix< T >.

Definition at line 492 of file petsc_matrix.C.

493{
495}

◆ l1_norm_diff()

template<typename T >
Real libMesh::SparseMatrix< T >::l1_norm_diff ( const SparseMatrix< T > &  other_mat) const
inherited
Returns
The l1_norm() of the difference of this and other_mat

Definition at line 1362 of file sparse_matrix.C.

1363{
1364 auto diff_mat = this->clone();
1365 diff_mat->add(-1.0, other_mat);
1366 return diff_mat->l1_norm();
1367}
virtual std::unique_ptr< SparseMatrix< T > > clone() const =0

Referenced by libMesh::l1_norm_diff().

◆ linfty_norm()

template<typename T >
Real libMesh::PetscMatrix< T >::linfty_norm ( ) const
overridevirtual
Returns
The \( \ell_{\infty} \)-norm of the matrix, that is the max row sum:

\( |M|_{\infty} = \max_{i} \sum_{j} |M_{ij}| \)

This is the natural matrix norm that is compatible to the \( \ell_{\infty} \)-norm of vectors, i.e. \( |Mv|_{\infty} \leq |M|_{\infty} |v|_{\infty} \). (cf. Haemmerlin-Hoffmann : Numerische Mathematik)

Implements libMesh::SparseMatrix< T >.

Definition at line 502 of file petsc_matrix.C.

503{
505}

◆ local_m()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::local_m ( ) const
finalvirtualinherited

Get the number of rows owned by this process.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 142 of file petsc_matrix_base.C.

143{
144 libmesh_assert (this->initialized());
145
146 PetscInt m = 0;
147
148 LibmeshPetscCall(MatGetLocalSize (this->_mat, &m, NULL));
149
150 return static_cast<numeric_index_type>(m);
151}

References libMesh::initialized(), and libMesh::libmesh_assert().

◆ local_n()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::local_n ( ) const
finalvirtualinherited

Get the number of columns owned by this process.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 166 of file petsc_matrix_base.C.

167{
168 libmesh_assert (this->initialized());
169
170 PetscInt n = 0;
171
172 LibmeshPetscCall(MatGetLocalSize (this->_mat, NULL, &n));
173
174 return static_cast<numeric_index_type>(n);
175}

◆ m()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::m ( ) const
overridevirtualinherited
Returns
The row-dimension of the matrix.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 130 of file petsc_matrix_base.C.

131{
132 libmesh_assert (this->initialized());
133
134 PetscInt petsc_m=0, petsc_n=0;
135
136 LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
137
138 return static_cast<numeric_index_type>(petsc_m);
139}

◆ mat() [1/2]

template<typename T >
Mat libMesh::PetscMatrixBase< T >::mat ( )
inlineinherited

◆ mat() [2/2]

template<typename T >
Mat libMesh::PetscMatrixBase< T >::mat ( ) const
inlineinherited

Definition at line 121 of file petsc_matrix_base.h.

121{ libmesh_assert(_mat); return _mat; }

References libMesh::PetscMatrixBase< T >::_mat, and libMesh::libmesh_assert().

◆ matrix_matrix_mult()

template<typename T >
void libMesh::PetscMatrix< T >::matrix_matrix_mult ( SparseMatrix< T > &  X,
SparseMatrix< T > &  Y,
bool  reuse = false 
)
overridevirtual

Compute Y = A*X for matrix X.

Set reuse = true if this->_mat and X have the same nonzero pattern as before, and Y is obtained from a previous call to this function with reuse = false

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1043 of file petsc_matrix.C.

1044{
1045 libmesh_assert (this->initialized());
1046
1047 // sanity check
1048 // we do not check the Y_out size here as we will initialize & close it at the end
1049 libmesh_assert_equal_to (this->n(), X_in.m());
1050
1051 const PetscMatrix<T> * X = cast_ptr<const PetscMatrix<T> *> (&X_in);
1052 PetscMatrix<T> * Y = cast_ptr<PetscMatrix<T> *> (&Y_out);
1053
1054 // the matrix from which we copy the values has to be assembled/closed
1055 libmesh_assert(X->closed());
1056
1057 semiparallel_only();
1058
1059 if (reuse)
1060 LibmeshPetscCall(MatMatMult(this->_mat, X->_mat, MAT_REUSE_MATRIX, PETSC_DEFAULT, &Y->_mat));
1061 else
1062 {
1063 Y->clear();
1064 LibmeshPetscCall(MatMatMult(this->_mat, X->_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Y->_mat));
1065 }
1066 // Specify that the new matrix is initialized
1067 // We do not close it here as `MatMatMult` ensures Y being closed
1068 Y->_is_initialized = true;
1069}

References libMesh::PetscMatrixBase< T >::_mat, libMesh::PetscMatrixBase< T >::closed(), libMesh::initialized(), libMesh::libmesh_assert(), and libMesh::SparseMatrix< T >::m().

Referenced by libMesh::PetscPreconditioner< T >::set_hypre_ads_data().

◆ n()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::n ( ) const
overridevirtualinherited
Returns
The column-dimension of the matrix.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 154 of file petsc_matrix_base.C.

155{
156 libmesh_assert (this->initialized());
157
158 PetscInt petsc_m=0, petsc_n=0;
159
160 LibmeshPetscCall(MatGetSize (this->_mat, &petsc_m, &petsc_n));
161
162 return static_cast<numeric_index_type>(petsc_n);
163}

◆ n_nonzeros()

template<typename T >
std::size_t libMesh::SparseMatrix< T >::n_nonzeros ( ) const
virtualinherited
Returns
the global number of non-zero entries in the matrix sparsity pattern

Definition at line 266 of file sparse_matrix.C.

267{
268 if (!_sp)
269 return 0;
270 return _sp->n_nonzeros();
271}
std::size_t n_nonzeros() const
The total number of nonzeros in the global matrix.

◆ n_objects()

static unsigned int libMesh::ReferenceCounter::n_objects ( )
inlinestaticinherited

Prints the number of outstanding (created, but not yet destroyed) objects.

Definition at line 85 of file reference_counter.h.

86 { return _n_objects; }
static Threads::atomic< unsigned int > _n_objects
The number of objects.

References libMesh::ReferenceCounter::_n_objects.

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ n_processors()

processor_id_type libMesh::ParallelObject::n_processors ( ) const
inlineinherited
Returns
The number of processors in the group.

Definition at line 103 of file parallel_object.h.

104 {
105 processor_id_type returnval =
106 cast_int<processor_id_type>(_communicator.size());
107 libmesh_assert(returnval); // We never have an empty comm
108 return returnval;
109 }
processor_id_type size() const
uint8_t processor_id_type
Definition id_types.h:104

References libMesh::ParallelObject::_communicator, libMesh::libmesh_assert(), and libMesh::Parallel::Communicator::size().

Referenced by libMesh::Partitioner::_find_global_index_by_pid_map(), libMesh::BoundaryInfo::_find_id_maps(), libMesh::DofMap::add_constraints_to_send_list(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::DistributedMesh::add_elem(), libMesh::DistributedMesh::add_node(), libMesh::System::add_vector(), libMesh::LaplaceMeshSmoother::allgather_graph(), libMesh::DofMap::allgather_recursive_constraints(), libMesh::FEMSystem::assembly(), libMesh::Nemesis_IO::assert_symmetric_cmaps(), libMesh::Partitioner::assign_partitioning(), libMesh::AztecLinearSolver< T >::AztecLinearSolver(), libMesh::Partitioner::build_graph(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::DistributedMesh::clear(), libMesh::DistributedMesh::clear_elems(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::UnstructuredMesh::copy_nodes_and_elements(), libMesh::ExodusII_IO::copy_scalar_solution(), libMesh::Nemesis_IO::copy_scalar_solution(), libMesh::UnstructuredMesh::create_pid_mesh(), libMesh::MeshTools::create_processor_bounding_box(), libMesh::DistributedMesh::DistributedMesh(), libMesh::EnsightIO::EnsightIO(), libMesh::RBEIMEvaluation::gather_bfs(), libMesh::MeshBase::get_info(), libMesh::StaticCondensation::init(), libMesh::SystemSubsetBySubdomain::init(), libMesh::PetscDMWrapper::init_petscdm(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::Nemesis_IO_Helper::initialize(), libMesh::DistributedMesh::insert_elem(), libMesh::NumericVector< T >::is_effectively_ghosted(), libMesh::NumericVector< T >::is_effectively_serial(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_new_node_procids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_topology_consistent_procids< Node >(), libMesh::MeshTools::libmesh_assert_valid_boundary_ids(), libMesh::MeshTools::libmesh_assert_valid_dof_ids(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::MeshTools::libmesh_assert_valid_refinement_flags(), libMesh::MeshBase::n_active_elem_on_proc(), libMesh::DofMap::n_dofs_per_processor(), libMesh::MeshBase::n_elem_on_proc(), libMesh::MeshBase::n_nodes_on_proc(), libMesh::RBEIMEvaluation::node_gather_bfs(), libMesh::MeshBase::partition(), libMesh::Partitioner::partition(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::Partitioner::partition_unpartitioned_elements(), libMesh::MeshBase::print_constraint_rows(), libMesh::DofMap::print_dof_constraints(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::NameBasedIO::read(), libMesh::CheckpointIO::read_connectivity(), libMesh::XdrIO::read_header(), libMesh::CheckpointIO::read_nodes(), libMesh::System::read_parallel_data(), libMesh::System::read_SCALAR_dofs(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::System::read_serialized_vector(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::Partitioner::repartition(), OverlappingFunctorTest::run_partitioner_test(), libMesh::DofMap::scatter_constraints(), libMesh::DistributedMesh::set_next_unique_id(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), WriteVecAndScalar::setupTests(), libMesh::RBEIMEvaluation::side_gather_bfs(), DistributedMeshTest::testRemoteElemError(), CheckpointIOTest::testSplitter(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::GMVIO::write_binary(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::VTKIO::write_nodal_data(), libMesh::System::write_parallel_data(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::XdrIO::write_serialized_nodes(), and libMesh::XdrIO::write_serialized_nodesets().

◆ need_full_sparsity_pattern()

template<typename T >
virtual bool libMesh::SparseMatrix< T >::need_full_sparsity_pattern ( ) const
inlinevirtualinherited
Returns
true if this sparse matrix format needs to be fed the graph of the sparse matrix.

This is true for LaspackMatrix, but not PetscMatrixBase subclasses. In the case where the full graph is not required, we can efficiently approximate it to provide a good estimate of the required size of the sparse matrix.

Reimplemented in libMesh::LaspackMatrix< T >, and libMesh::EpetraMatrix< T >.

Definition at line 162 of file sparse_matrix.h.

163 { return false; }

Referenced by libMesh::DofMap::attach_matrix(), and libMesh::DofMap::update_sparsity_pattern().

◆ norm()

template<typename T >
template<NormType N>
Real libMesh::PetscMatrix< T >::norm ( ) const
private
Returns
A norm of the matrix, where the type of norm to compute is determined by the template parameter N of the PETSc-defined type NormType. The valid template arguments are NORM_1, NORM_FROBENIUS and NORM_INFINITY, as used to define l1_norm(), frobenius_norm() and linfty_norm().

Definition at line 474 of file petsc_matrix.C.

475{
476 libmesh_assert (this->initialized());
477
478 semiparallel_only();
479
480 PetscReal petsc_value;
481 Real value;
482
483 libmesh_assert (this->closed());
484
485 LibmeshPetscCall(MatNorm(this->_mat, N, &petsc_value));
486
487 value = static_cast<Real>(petsc_value);
488
489 return value;
490}
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real

References libMesh::closed(), libMesh::initialized(), libMesh::libmesh_assert(), libMesh::Real, and value.

◆ operator()()

template<typename T >
T libMesh::PetscMatrix< T >::operator() ( const numeric_index_type  i,
const numeric_index_type  j 
) const
overridevirtual
Returns
A copy of matrix entry (i,j).
Note
This may be an expensive operation, and you should always be careful where you call this function.

Implements libMesh::SparseMatrix< T >.

Definition at line 1123 of file petsc_matrix.C.

1125{
1126 libmesh_assert (this->initialized());
1127
1128 // If the entry is not in the sparse matrix, it is 0.
1129 T value=0.;
1130
1131 PetscInt
1132 i_val=static_cast<PetscInt>(i_in),
1133 j_val=static_cast<PetscInt>(j_in);
1134
1135
1136 // the matrix needs to be closed for this to work
1137 // this->close();
1138 // but closing it is a semiparallel operation; we want operator()
1139 // to run on one processor.
1140 libmesh_assert(this->closed());
1141
1142 LibmeshPetscCall(MatGetValue(this->_mat, i_val, j_val, &value));
1143
1144 return value;
1145}

References libMesh::closed(), libMesh::initialized(), libMesh::libmesh_assert(), and value.

◆ operator=() [1/3]

template<typename T >
PetscMatrix< T > & libMesh::PetscMatrix< T >::operator= ( const PetscMatrix< T > &  v)

Definition at line 1198 of file petsc_matrix.C.

1199{
1200 semiparallel_only();
1201
1202 if (this->_mat)
1203 {
1204 PetscBool assembled;
1205 LibmeshPetscCall(MatAssembled(this->_mat, &assembled));
1206#ifndef NDEBUG
1207 const bool cxx_assembled = (assembled == PETSC_TRUE) ? true : false;
1208 libmesh_assert(this->_communicator.verify(cxx_assembled));
1209#endif
1210
1211 if (!assembled)
1212 // MatCopy does not work with an unassembled matrix. We could use MatDuplicate but then we
1213 // would have to destroy the matrix we manage and others might be relying on that data. So
1214 // we just assemble here regardless of the preceding level of matrix fill
1215 this->close();
1216 LibmeshPetscCall(MatCopy(v._mat, this->_mat, DIFFERENT_NONZERO_PATTERN));
1217 }
1218 else
1219 LibmeshPetscCall(MatDuplicate(v._mat, MAT_COPY_VALUES, &this->_mat));
1220
1221 this->_is_initialized = true;
1222
1223 return *this;
1224}
timpi_pure bool verify(const T &r) const

References libMesh::PetscMatrixBase< T >::_mat, and libMesh::libmesh_assert().

◆ operator=() [2/3]

template<typename T >
SparseMatrix< T > & libMesh::PetscMatrix< T >::operator= ( const SparseMatrix< T > &  )
overridevirtual

This looks like a copy assignment operator, but note that, unlike normal copy assignment operators, it is pure virtual.

This function should be overridden in derived classes so that they can be copied correctly via references to the base class. This design usually isn't a good idea in general, but in this context it works because we usually don't have a mix of different kinds of SparseMatrix active in the library at a single time.

Returns
A reference to *this as the base type.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1227 of file petsc_matrix.C.

1228{
1229 *this = cast_ref<const PetscMatrix<T> &>(v);
1230 return *this;
1231}

◆ operator=() [3/3]

template<typename T >
PetscMatrix & libMesh::PetscMatrix< T >::operator= ( PetscMatrix< T > &&  )
delete

◆ preallocate()

template<typename T >
void libMesh::PetscMatrix< T >::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 
)
protected

Definition at line 260 of file petsc_matrix.C.

264{
265 // Make sure the sparsity pattern isn't empty unless the matrix is 0x0
266 libmesh_assert_equal_to (n_nz.size(), m_l);
267 libmesh_assert_equal_to (n_oz.size(), m_l);
268 // Avoid unused warnings when not configured with block storage
269 libmesh_ignore(blocksize_in);
270
271#ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
272 PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
273
274 if (blocksize > 1)
275 {
276 // transform the per-entry n_nz and n_oz arrays into their block counterparts.
277 std::vector<numeric_index_type> b_n_nz, b_n_oz;
278
279 transform_preallocation_arrays (blocksize,
280 n_nz, n_oz,
281 b_n_nz, b_n_oz);
282
283 LibmeshPetscCall(MatSeqBAIJSetPreallocation (this->_mat,
284 blocksize,
285 0,
286 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data())));
287
288 LibmeshPetscCall(MatMPIBAIJSetPreallocation (this->_mat,
289 blocksize,
290 0,
291 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()),
292 0,
293 numeric_petsc_cast(b_n_oz.empty() ? nullptr : b_n_oz.data())));
294 }
295 else
296#endif
297 {
298 switch (this->_mat_type) {
299 case AIJ:
300 LibmeshPetscCall(MatSeqAIJSetPreallocation (this->_mat,
301 0,
302 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data())));
303 LibmeshPetscCall(MatMPIAIJSetPreallocation (this->_mat,
304 0,
305 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
306 0,
307 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data())));
308 break;
309
310 case HYPRE:
311#if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
312 LibmeshPetscCall(MatHYPRESetPreallocation (this->_mat,
313 0,
314 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
315 0,
316 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data())));
317#else
318 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
319#endif
320 break;
321
322 default: libmesh_error_msg("Unsupported petsc matrix type");
323 }
324
325 }
326}

References libMesh::AIJ, libMesh::HYPRE, libMesh::libmesh_ignore(), and libMesh::numeric_petsc_cast().

◆ print() [1/3]

template<typename T >
void libMesh::SparseMatrix< T >::print ( const std::string &  filename) const
inherited

Print the contents of the matrix to a file, with a file format depending on the extension of filename.

Definition at line 754 of file sparse_matrix.C.

755{
756 {
757 std::ofstream outstr (filename.c_str());
758 libmesh_error_msg_if
759 (!outstr.good(), "ERROR: cannot write to file:\n\t" <<
760 filename);
761 }
762
763 std::string_view basename = Utility::basename_of(filename);
764
765 const bool gzipped_file = Utility::ends_with(filename, ".gz");
766
767 if (gzipped_file)
768 basename.remove_suffix(3);
769
770 if (Utility::ends_with(basename, ".matlab") ||
771 Utility::ends_with(basename, ".m"))
772 this->print_matlab(filename);
773 else if (Utility::ends_with(basename, ".petsc64"))
774 {
775#ifndef LIBMESH_HAVE_PETSC
776 libmesh_error_msg("Cannot load PETSc matrix file " <<
777 filename << " without PETSc-enabled libMesh.");
778#elif LIBMESH_DOF_ID_BYTES != 8
779 libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
780 filename << " with non-64-bit libMesh.");
781#endif
782 if (gzipped_file)
783 libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
784 this->print_petsc_binary(filename);
785 }
786 else if (Utility::ends_with(basename, ".petsc32"))
787 {
788#ifndef LIBMESH_HAVE_PETSC
789 libmesh_error_msg("Cannot load PETSc matrix file " <<
790 filename << " without PETSc-enabled libMesh.");
791#elif LIBMESH_DOF_ID_BYTES != 4
792 libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
793 filename << " with non-32-bit libMesh.");
794#endif
795 if (gzipped_file)
796 libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
797 this->print_petsc_binary(filename);
798 }
799 else if (Utility::ends_with(basename, ".h5"))
800 {
801 if (gzipped_file)
802 libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
803 this->print_coreform_hdf5(filename);
804 }
805 else
806 libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
807 << basename
808 << "\n I understand the following:\n\n"
809 << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
810 << " *.matlab -- Matlab sparse matrix format\n"
811 << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
812 << " *.m -- Matlab sparse matrix format\n"
813 << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
814 << " *.petsc32 -- PETSc binary format, 32-bit\n"
815 << " *.petsc64 -- PETSc binary format, 64-bit\n"
816 );
817}
virtual void print_coreform_hdf5(const std::string &filename, const std::string &groupname="extraction") const
Print the contents of the matrix to a file, with the HDF5 sparse matrix format used by CoreForm,...
virtual void print_petsc_binary(const std::string &filename) const
Write the contents of the matrix to a file in PETSc's binary sparse matrix format.
virtual void print_matlab(const std::string &="") const
Print the contents of the matrix in Matlab's sparse matrix format.
bool ends_with(std::string_view superstring, std::string_view suffix)
Look for a substring at the very end of a string.
Definition utility.C:213
std::string_view basename_of(const std::string &fullname)
Definition utility.C:108

References libMesh::Utility::basename_of(), and libMesh::Utility::ends_with().

◆ print() [2/3]

void libMesh::SparseMatrix< Complex >::print ( std::ostream &  os,
const bool  sparse 
) const
inherited

Definition at line 159 of file sparse_matrix.C.

160{
161 // std::complex<>::operator<<() is defined, but use this form
162
163 if (sparse)
164 {
165 libmesh_not_implemented();
166 }
167
168 os << "Real part:" << std::endl;
169 for (auto i : make_range(this->m()))
170 {
171 for (auto j : make_range(this->n()))
172 os << std::setw(8) << (*this)(i,j).real() << " ";
173 os << std::endl;
174 }
175
176 os << std::endl << "Imaginary part:" << std::endl;
177 for (auto i : make_range(this->m()))
178 {
179 for (auto j : make_range(this->n()))
180 os << std::setw(8) << (*this)(i,j).imag() << " ";
181 os << std::endl;
182 }
183}
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition int_range.h:176
boost::multiprecision::float128 real(const boost::multiprecision::float128 in)
boost::multiprecision::float128 imag(const boost::multiprecision::float128)

References libMesh::make_range().

◆ print() [3/3]

template<typename T >
void libMesh::SparseMatrix< T >::print ( std::ostream &  os = libMesh::out,
const bool  sparse = false 
) const
inherited

Print the contents of the matrix to the screen in a uniform style, regardless of matrix/solver package being used.

Definition at line 275 of file sparse_matrix.C.

276{
277 parallel_object_only();
278
279 libmesh_assert (this->initialized());
280
281 const numeric_index_type first_dof = this->row_start(),
282 end_dof = this->row_stop();
283
284 // We'll print the matrix from processor 0 to make sure
285 // it's serialized properly
286 if (this->processor_id() == 0)
287 {
288 libmesh_assert_equal_to (first_dof, 0);
289 for (numeric_index_type i : make_range(end_dof))
290 {
291 if (sparse)
292 {
293 for (auto j : make_range(this->n()))
294 {
295 T c = (*this)(i,j);
296 if (c != static_cast<T>(0.0))
297 {
298 os << i << " " << j << " " << c << std::endl;
299 }
300 }
301 }
302 else
303 {
304 for (auto j : make_range(this->n()))
305 os << (*this)(i,j) << " ";
306 os << std::endl;
307 }
308 }
309
310 std::vector<numeric_index_type> ibuf, jbuf;
311 std::vector<T> cbuf;
312 numeric_index_type currenti = end_dof;
313 for (auto p : IntRange<processor_id_type>(1, this->n_processors()))
314 {
315 this->comm().receive(p, ibuf);
316 this->comm().receive(p, jbuf);
317 this->comm().receive(p, cbuf);
318 libmesh_assert_equal_to (ibuf.size(), jbuf.size());
319 libmesh_assert_equal_to (ibuf.size(), cbuf.size());
320
321 if (ibuf.empty())
322 continue;
323 libmesh_assert_greater_equal (ibuf.front(), currenti);
324 libmesh_assert_greater_equal (ibuf.back(), ibuf.front());
325
326 std::size_t currentb = 0;
327 for (;currenti <= ibuf.back(); ++currenti)
328 {
329 if (sparse)
330 {
331 for (numeric_index_type j=0; j<this->n(); j++)
332 {
333 if (currentb < ibuf.size() &&
334 ibuf[currentb] == currenti &&
335 jbuf[currentb] == j)
336 {
337 os << currenti << " " << j << " " << cbuf[currentb] << std::endl;
338 currentb++;
339 }
340 }
341 }
342 else
343 {
344 for (auto j : make_range(this->n()))
345 {
346 if (currentb < ibuf.size() &&
347 ibuf[currentb] == currenti &&
348 jbuf[currentb] == j)
349 {
350 os << cbuf[currentb] << " ";
351 currentb++;
352 }
353 else
354 os << static_cast<T>(0.0) << " ";
355 }
356 os << std::endl;
357 }
358 }
359 }
360 if (!sparse)
361 {
362 for (; currenti != this->m(); ++currenti)
363 {
364 for (numeric_index_type j=0; j<this->n(); j++)
365 os << static_cast<T>(0.0) << " ";
366 os << std::endl;
367 }
368 }
369 }
370 else
371 {
372 std::vector<numeric_index_type> ibuf, jbuf;
373 std::vector<T> cbuf;
374
375 // We'll assume each processor has access to entire
376 // matrix rows, so (*this)(i,j) is valid if i is a local index.
377 for (numeric_index_type i : make_range(first_dof, end_dof))
378 {
379 for (auto j : make_range(this->n()))
380 {
381 T c = (*this)(i,j);
382 if (c != static_cast<T>(0.0))
383 {
384 ibuf.push_back(i);
385 jbuf.push_back(j);
386 cbuf.push_back(c);
387 }
388 }
389 }
390 this->comm().send(0,ibuf);
391 this->comm().send(0,jbuf);
392 this->comm().send(0,cbuf);
393 }
394}
Status receive(const unsigned int dest_processor_id, T &buf, const MessageTag &tag=any_tag) const
void send(const unsigned int dest_processor_id, const T &buf, const MessageTag &tag=no_tag) const
processor_id_type processor_id() const
processor_id_type n_processors() const
virtual numeric_index_type row_stop() const =0
virtual numeric_index_type row_start() const =0

References libMesh::initialized(), libMesh::libmesh_assert(), and libMesh::make_range().

Referenced by libMesh::EigenSparseMatrix< T >::print_personal(), and libMesh::LaspackMatrix< T >::print_personal().

◆ print_coreform_hdf5()

template<typename T >
void libMesh::SparseMatrix< T >::print_coreform_hdf5 ( const std::string &  filename,
const std::string &  groupname = "extraction" 
) const
virtualinherited

Print the contents of the matrix to a file, with the HDF5 sparse matrix format used by CoreForm, putting CSR sparse matrix data in the group given by groupname.

Definition at line 521 of file sparse_matrix.C.

523{
524#if defined(LIBMESH_USE_COMPLEX_NUMBERS) || !defined(LIBMESH_HAVE_HDF5)
525 // TODO: HDF5 version 2.0.0 and later adds native support for
526 // complex numbers with the H5T_NATIVE_DOUBLE_COMPLEX type [0], so
527 // we could consider supporting T==std::complex<Real> in the future.
528 // [0]: https://forum.hdfgroup.org/t/coming-in-the-next-hdf5-release-native-support-for-complex-number-datatypes/13543
529 libmesh_ignore(filename, groupname);
530 libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
531#else
532 LOG_SCOPE("print_coreform_hdf5()", "SparseMatrix");
533
534 // In this implementation, we copy the SparseMatrix entries into a
535 // std::vector<double>, so this won't work for any Number type for
536 // which sizeof(Number) > sizeof(double).
537 if constexpr (sizeof(T) > sizeof(double))
538 libmesh_not_implemented();
539
540 const numeric_index_type first_dof = this->row_start(),
541 end_dof = this->row_stop();
542
543 std::vector<std::size_t> cols, row_offsets;
544 std::vector<double> vals;
545
546 if (this->processor_id() == 0)
547 row_offsets.push_back(0);
548
549 for (numeric_index_type i : make_range(first_dof, end_dof))
550 {
551 for (auto j : make_range(this->n()))
552 {
553 T c = (*this)(i,j);
554 if (c != static_cast<T>(0.0))
555 {
556 cols.push_back(j);
557 vals.push_back(c);
558 }
559 }
560 // This is a *local* row offset; proc 0 may need to adjust later
561 row_offsets.push_back(cols.size());
562 }
563
564 if (this->processor_id() == 0)
565 {
566 const hid_t file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC,
567 H5P_DEFAULT, H5P_DEFAULT);
568 if (file == H5I_INVALID_HID)
569 libmesh_file_error(filename);
570
571 const hid_t group =
572 H5Gcreate2(file, groupname.c_str(), H5P_DEFAULT, H5P_DEFAULT,
573 H5P_DEFAULT);
574 check_open(filename, group, groupname);
575
576 auto write_size_attribute = [&filename, &group]
577 (const std::string & attribute_name, unsigned long long writeval)
578 {
579 const hid_t fspace = H5Screate(H5S_SCALAR);
580 check_hdf5(filename, fspace, attribute_name + " fspace");
581
582 const hid_t attr = H5Acreate2(group, attribute_name.c_str(),
583 H5T_STD_I64LE, fspace,
584 H5P_DEFAULT, H5P_DEFAULT);
585 check_hdf5(filename, attr, attribute_name);
586
587 // HDF5 is supposed to handle both upscaling and endianness
588 // conversions here
589 const herr_t errval = H5Awrite(attr, H5T_NATIVE_ULLONG, &writeval);
590 check_hdf5(filename, errval, attribute_name + " write");
591
592 H5Aclose(attr);
593 H5Sclose(fspace);
594 };
595
596 write_size_attribute("num_cols", this->n());
597 write_size_attribute("num_rows", this->m());
598
599 auto write_vector = [&filename, &group]
600 (const std::string & dataname, auto hdf5_file_type,
601 auto hdf5_native_type, auto & datavec)
602 {
603 const hsize_t len[1] = {cast_int<hsize_t>(datavec.size())};
604
605 const hid_t space = H5Screate_simple(1, len, nullptr);
606 check_hdf5(filename, space, dataname + " space");
607
608 const hid_t data =
609 H5Dcreate2(group, dataname.c_str(), hdf5_file_type, space,
610 H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
611 check_hdf5(filename, data, dataname + " data");
612
613 const hid_t errval =
614 H5Dwrite(data, hdf5_native_type, H5S_ALL, H5S_ALL,
615 H5P_DEFAULT, datavec.data());
616 check_hdf5(filename, errval, dataname + " write");
617
618 H5Dclose(data);
619 H5Sclose(space);
620 };
621
622 std::vector<std::size_t> vals_sizes, first_dofs;
623 this->comm().gather(0, vals.size(), vals_sizes);
624 this->comm().gather(0, cast_int<std::size_t>(first_dof), first_dofs);
625 first_dofs.push_back(this->m());
626
627 this->comm().allgather(cols);
628 this->comm().allgather(vals);
629 this->comm().allgather(row_offsets);
630
631 libmesh_assert_equal_to(vals.size(),
632 cols.size());
633 libmesh_assert_equal_to(vals.size(),
634 std::accumulate(vals_sizes.begin(),
635 vals_sizes.end(),
636 std::size_t(0)));
637
638 std::size_t extra_offset = 0;
639 for (auto p : make_range(processor_id_type(1), this->n_processors()))
640 {
641 extra_offset += vals_sizes[p-1];
642 for (auto i : make_range(first_dofs[p]+1, first_dofs[p+1]+1))
643 row_offsets[i] += extra_offset;
644 }
645
646 write_vector("cols", H5T_STD_U64LE, H5T_NATIVE_ULLONG, cols);
647 write_vector("row_offsets", H5T_STD_U64LE, H5T_NATIVE_ULLONG, row_offsets);
648 write_vector("vals", H5T_IEEE_F64LE, H5T_NATIVE_DOUBLE, vals);
649
650 H5Gclose(group);
651 H5Fclose(file);
652 }
653 else
654 {
655 std::vector<std::size_t> dummy;
656 this->comm().gather(0, vals.size(), dummy);
657 this->comm().gather(0, cast_int<std::size_t>(first_dof), dummy);
658 this->comm().allgather(cols);
659 this->comm().allgather(vals);
660 this->comm().allgather(row_offsets);
661 }
662#endif // LIBMESH_HAVE_HDF5
663}
void gather(const unsigned int root_id, const T &send_data, std::vector< T, A > &recv) const
void allgather(const T &send_data, std::vector< T, A > &recv_data) const

◆ print_info()

void libMesh::ReferenceCounter::print_info ( std::ostream &  out_stream = libMesh::out)
staticinherited

Prints the reference information, by default to libMesh::out.

Definition at line 81 of file reference_counter.C.

82{
84 out_stream << ReferenceCounter::get_info();
85}
static std::string get_info()
Gets a string containing the reference information.

References libMesh::ReferenceCounter::_enable_print_counter, and libMesh::ReferenceCounter::get_info().

Referenced by libMesh::LibMeshInit::~LibMeshInit().

◆ print_matlab()

template<typename T >
void libMesh::PetscMatrix< T >::print_matlab ( const std::string &  name = "") const
overridevirtual

Print the contents of the matrix in Matlab's sparse matrix format.

Optionally prints the matrix to the file named name. If name is not specified it is dumped to the screen.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 510 of file petsc_matrix.C.

511{
512 libmesh_assert (this->initialized());
513
514 semiparallel_only();
515
516 if (!this->closed())
517 {
518 libmesh_deprecated();
519 libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_matlab().\n"
520 "Please update your code, as this warning will become an error in a future release.");
521 const_cast<PetscMatrix<T> *>(this)->close();
522 }
523
524 // Create an ASCII file containing the matrix
525 // if a filename was provided.
526 if (name != "")
527 {
528 WrappedPetsc<PetscViewer> petsc_viewer;
529
530 LibmeshPetscCall(PetscViewerASCIIOpen( this->comm().get(),
531 name.c_str(),
532 petsc_viewer.get()));
533
534#if PETSC_VERSION_LESS_THAN(3,7,0)
535 LibmeshPetscCall(PetscViewerSetFormat (petsc_viewer,
536 PETSC_VIEWER_ASCII_MATLAB));
537#else
538 LibmeshPetscCall(PetscViewerPushFormat (petsc_viewer,
539 PETSC_VIEWER_ASCII_MATLAB));
540#endif
541
542 LibmeshPetscCall(MatView (this->_mat, petsc_viewer));
543 }
544
545 // Otherwise the matrix will be dumped to the screen.
546 else
547 {
548#if PETSC_VERSION_LESS_THAN(3,7,0)
549 LibmeshPetscCall(PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
550 PETSC_VIEWER_ASCII_MATLAB));
551#else
552 LibmeshPetscCall(PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
553 PETSC_VIEWER_ASCII_MATLAB));
554#endif
555
556 LibmeshPetscCall(MatView (this->_mat, PETSC_VIEWER_STDOUT_WORLD));
557 }
558}

References libMesh::closed(), libMesh::WrappedPetsc< T >::get(), libMesh::initialized(), and libMesh::libmesh_assert().

◆ print_personal()

template<typename T >
void libMesh::PetscMatrix< T >::print_personal ( std::ostream &  os = libMesh::out) const
overridevirtual

Print the contents of the matrix to the screen with the PETSc viewer.

This function only allows printing to standard out since we have limited ourselves to one PETSc implementation for writing.

Implements libMesh::SparseMatrix< T >.

Definition at line 565 of file petsc_matrix.C.

566{
567 libmesh_assert (this->initialized());
568
569 // Routine must be called in parallel on parallel matrices
570 // and serial on serial matrices.
571 semiparallel_only();
572
573 // #ifndef NDEBUG
574 // if (os != std::cout)
575 // libMesh::err << "Warning! PETSc can only print to std::cout!" << std::endl;
576 // #endif
577
578 // Matrix must be in an assembled state to be printed
579 if (!this->closed())
580 {
581 libmesh_deprecated();
582 libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_personal().\n"
583 "Please update your code, as this warning will become an error in a future release.");
584 const_cast<PetscMatrix<T> *>(this)->close();
585 }
586
587 // Print to screen if ostream is stdout
588 if (os.rdbuf() == std::cout.rdbuf())
589 LibmeshPetscCall(MatView(this->_mat, NULL));
590
591 // Otherwise, print to the requested file, in a roundabout way...
592 else
593 {
594 // We will create a temporary filename, and file, for PETSc to
595 // write to.
596 std::string temp_filename;
597
598 {
599 // Template for temporary filename
600 char c[] = "temp_petsc_matrix.XXXXXX";
601
602 // Generate temporary, unique filename only on processor 0. We will
603 // use this filename for PetscViewerASCIIOpen, before copying it into
604 // the user's stream
605 if (this->processor_id() == 0)
606 {
607 int fd = mkstemp(c);
608
609 // Check to see that mkstemp did not fail.
610 libmesh_error_msg_if(fd == -1, "mkstemp failed in PetscMatrix::print_personal()");
611
612 // mkstemp returns a file descriptor for an open file,
613 // so let's close it before we hand it to PETSc!
614 ::close (fd);
615 }
616
617 // Store temporary filename as string, makes it easier to broadcast
618 temp_filename = c;
619 }
620
621 // Now broadcast the filename from processor 0 to all processors.
622 this->comm().broadcast(temp_filename);
623
624 // PetscViewer object for passing to MatView
625 PetscViewer petsc_viewer;
626
627 // This PETSc function only takes a string and handles the opening/closing
628 // of the file internally. Since print_personal() takes a reference to
629 // an ostream, we have to do an extra step... print_personal() should probably
630 // have a version that takes a string to get rid of this problem.
631 LibmeshPetscCall(PetscViewerASCIIOpen( this->comm().get(),
632 temp_filename.c_str(),
633 &petsc_viewer));
634
635 // Probably don't need to set the format if it's default...
636 // ierr = PetscViewerSetFormat (petsc_viewer,
637 // PETSC_VIEWER_DEFAULT);
638 // LIBMESH_CHKERR(ierr);
639
640 // Finally print the matrix using the viewer
641 LibmeshPetscCall(MatView (this->_mat, petsc_viewer));
642
643 if (this->processor_id() == 0)
644 {
645 // Now the inefficient bit: open temp_filename as an ostream and copy the contents
646 // into the user's desired ostream. We can't just do a direct file copy, we don't even have the filename!
647 std::ifstream input_stream(temp_filename.c_str());
648 os << input_stream.rdbuf(); // The "most elegant" way to copy one stream into another.
649 // os.close(); // close not defined in ostream
650
651 // Now remove the temporary file
652 input_stream.close();
653 std::remove(temp_filename.c_str());
654 }
655 }
656}
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
int mkstemp(char *tmpl)
Definition win_mkstemp.h:13

References libMesh::closed(), libMesh::initialized(), libMesh::libmesh_assert(), and mkstemp().

◆ print_petsc_binary()

template<typename T >
void libMesh::PetscMatrix< T >::print_petsc_binary ( const std::string &  filename) const
overridevirtual

Write the contents of the matrix to a file in PETSc's binary sparse matrix format.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 691 of file petsc_matrix.C.

692{
693 libmesh_assert (this->initialized());
694
695 // Our helper here isn't const-correct for writes
696 PetscMatrix<T> * nonconst_this = const_cast<PetscMatrix<T>*>(this);
697 nonconst_this->_petsc_viewer(filename, PETSCVIEWERBINARY, FILE_MODE_WRITE);
698}
void _petsc_viewer(const std::string &filename, PetscViewerType viewertype, PetscFileMode filemode)

References libMesh::PetscMatrix< T >::_petsc_viewer(), libMesh::initialized(), and libMesh::libmesh_assert().

◆ print_petsc_hdf5()

template<typename T >
void libMesh::PetscMatrix< T >::print_petsc_hdf5 ( const std::string &  filename) const
overridevirtual

Write the contents of the matrix to a file in PETSc's HDF5 sparse matrix format.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 703 of file petsc_matrix.C.

704{
705 libmesh_assert (this->initialized());
706
707 // Our helper here isn't const-correct for writes
708 PetscMatrix<T> * nonconst_this = const_cast<PetscMatrix<T>*>(this);
709 nonconst_this->_petsc_viewer(filename, PETSCVIEWERHDF5, FILE_MODE_WRITE);
710}

References libMesh::PetscMatrix< T >::_petsc_viewer(), libMesh::initialized(), and libMesh::libmesh_assert().

◆ processor_id()

processor_id_type libMesh::ParallelObject::processor_id ( ) const
inlineinherited
Returns
The rank of this processor in the group.

Definition at line 114 of file parallel_object.h.

115 { return cast_int<processor_id_type>(_communicator.rank()); }
processor_id_type rank() const

References libMesh::ParallelObject::_communicator, and libMesh::Parallel::Communicator::rank().

Referenced by libMesh::BoundaryInfo::_find_id_maps(), libMesh::PetscDMWrapper::add_dofs_to_section(), libMesh::DistributedMesh::add_elem(), libMesh::BoundaryInfo::add_elements(), libMesh::DistributedMesh::add_node(), libMesh::MeshTools::Modification::all_tri(), libMesh::FEMSystem::assembly(), libMesh::Nemesis_IO::assert_symmetric_cmaps(), libMesh::Partitioner::assign_partitioning(), libMesh::Nemesis_IO_Helper::build_element_and_node_maps(), libMesh::Partitioner::build_graph(), libMesh::InfElemBuilder::build_inf_elem(), libMesh::BoundaryInfo::build_node_list_from_side_list(), libMesh::EquationSystems::build_parallel_elemental_solution_vector(), libMesh::EquationSystems::build_parallel_solution_vector(), libMesh::MeshFunction::check_found_elem(), libMesh::DistributedMesh::clear(), libMesh::DistributedMesh::clear_elems(), libMesh::ExodusII_IO_Helper::close(), libMesh::Nemesis_IO_Helper::compute_border_node_ids(), libMesh::Nemesis_IO_Helper::compute_communication_map_parameters(), libMesh::Nemesis_IO_Helper::compute_internal_and_border_elems_and_internal_nodes(), libMesh::RBConstruction::compute_max_error_bound(), libMesh::Nemesis_IO_Helper::compute_node_communication_maps(), libMesh::Nemesis_IO_Helper::compute_num_global_elem_blocks(), libMesh::Nemesis_IO_Helper::compute_num_global_nodesets(), libMesh::Nemesis_IO_Helper::compute_num_global_sidesets(), libMesh::Nemesis_IO_Helper::construct_nemesis_filename(), libMesh::ExodusII_IO::copy_elemental_solution(), libMesh::ExodusII_IO::copy_nodal_solution(), libMesh::ExodusII_IO::copy_scalar_solution(), libMesh::Nemesis_IO::copy_scalar_solution(), libMesh::MeshTools::correct_node_proc_ids(), libMesh::ExodusII_IO_Helper::create(), libMesh::MeshCommunication::delete_remote_elements(), libMesh::DistributedMesh::DistributedMesh(), libMesh::DistributedMesh::DistributedMesh(), libMesh::DofMapBase::end_dof(), libMesh::DofMapBase::end_old_dof(), libMesh::EnsightIO::EnsightIO(), libMesh::GenericProjector< FFunctor, GFunctor, FValue, ProjectionAction >::SubFunctor::find_dofs_to_send(), libMesh::UnstructuredMesh::find_neighbors(), libMesh::DofMapBase::first_dof(), libMesh::DofMapBase::first_old_dof(), libMesh::RBEIMEvaluation::gather_bfs(), libMesh::Nemesis_IO_Helper::get_cmap_params(), libMesh::Nemesis_IO_Helper::get_eb_info_global(), libMesh::Nemesis_IO_Helper::get_elem_cmap(), libMesh::Nemesis_IO_Helper::get_elem_map(), libMesh::MeshBase::get_info(), libMesh::Nemesis_IO_Helper::get_init_global(), libMesh::Nemesis_IO_Helper::get_init_info(), libMesh::RBEIMEvaluation::get_interior_basis_functions_as_vecs(), libMesh::Nemesis_IO_Helper::get_loadbal_param(), libMesh::DofMap::get_local_constraints(), libMesh::MeshBase::get_local_constraints(), libMesh::Nemesis_IO_Helper::get_node_cmap(), libMesh::Nemesis_IO_Helper::get_node_map(), libMesh::Nemesis_IO_Helper::get_ns_param_global(), libMesh::Nemesis_IO_Helper::get_ss_param_global(), libMesh::SparsityPattern::Build::handle_vi_vj(), libMesh::LaplaceMeshSmoother::init(), libMesh::SystemSubsetBySubdomain::init(), HeatSystem::init_data(), libMesh::ExodusII_IO_Helper::initialize(), libMesh::ExodusII_IO_Helper::initialize_element_variables(), libMesh::ExodusII_IO_Helper::initialize_global_variables(), libMesh::ExodusII_IO_Helper::initialize_nodal_variables(), libMesh::DistributedMesh::insert_elem(), libMesh::MeshTools::Modification::interpolate_surface(), libMesh::SparsityPattern::Build::join(), libMesh::RBEvaluation::legacy_write_offline_data_to_files(), libMesh::RBSCMEvaluation::legacy_write_offline_data_to_files(), libMesh::TransientRBEvaluation::legacy_write_offline_data_to_files(), libMesh::MeshTools::libmesh_assert_consistent_distributed(), libMesh::MeshTools::libmesh_assert_consistent_distributed_nodes(), libMesh::MeshTools::libmesh_assert_contiguous_dof_ids(), libMesh::MeshTools::libmesh_assert_parallel_consistent_procids< Elem >(), libMesh::MeshTools::libmesh_assert_valid_neighbors(), libMesh::DistributedMesh::libmesh_assert_valid_parallel_object_ids(), main(), AugmentSparsityOnInterface::mesh_reinit(), libMesh::TriangulatorInterface::MeshedHole::MeshedHole(), libMesh::MeshBase::n_active_local_elem(), libMesh::BoundaryInfo::n_boundary_conds(), libMesh::MeshTools::n_connected_components(), libMesh::MeshBase::n_constraint_rows(), libMesh::BoundaryInfo::n_edge_conds(), libMesh::DofMapBase::n_local_dofs(), libMesh::MeshBase::n_local_elem(), libMesh::MeshBase::n_local_nodes(), libMesh::BoundaryInfo::n_nodeset_conds(), libMesh::BoundaryInfo::n_shellface_conds(), libMesh::RBEIMEvaluation::node_gather_bfs(), libMesh::DistributedMesh::own_node(), libMesh::BoundaryInfo::parallel_sync_node_ids(), libMesh::BoundaryInfo::parallel_sync_side_ids(), libMesh::MeshBase::print_constraint_rows(), libMesh::DofMap::print_dof_constraints(), libMesh::DofMap::process_mesh_constraint_rows(), libMesh::Nemesis_IO_Helper::put_cmap_params(), libMesh::Nemesis_IO_Helper::put_elem_cmap(), libMesh::Nemesis_IO_Helper::put_elem_map(), libMesh::Nemesis_IO_Helper::put_loadbal_param(), libMesh::Nemesis_IO_Helper::put_node_cmap(), libMesh::Nemesis_IO_Helper::put_node_map(), libMesh::XdrIO::read(), libMesh::Nemesis_IO::read(), libMesh::CheckpointIO::read(), libMesh::NameBasedIO::read(), libMesh::EquationSystems::read(), libMesh::EquationSystems::read(), libMesh::ExodusII_IO_Helper::read_elem_num_map(), libMesh::ExodusII_IO_Helper::read_global_values(), libMesh::CheckpointIO::read_header(), libMesh::ExodusII_IO::read_header(), libMesh::System::read_header(), libMesh::XdrIO::read_header(), libMesh::DynaIO::read_mesh(), libMesh::ExodusII_IO_Helper::read_node_num_map(), libMesh::System::read_parallel_data(), libMesh::RBConstruction::read_riesz_representors_from_files(), libMesh::TransientRBConstruction::read_riesz_representors_from_files(), libMesh::System::read_SCALAR_dofs(), libMesh::XdrIO::read_serialized_bc_names(), libMesh::XdrIO::read_serialized_bcs_helper(), libMesh::System::read_serialized_blocked_dof_objects(), libMesh::XdrIO::read_serialized_connectivity(), libMesh::System::read_serialized_data(), libMesh::XdrIO::read_serialized_nodes(), libMesh::XdrIO::read_serialized_nodesets(), libMesh::XdrIO::read_serialized_subdomain_names(), libMesh::System::read_serialized_vector(), libMesh::System::read_serialized_vectors(), libMesh::Nemesis_IO_Helper::read_var_names_impl(), libMesh::SimplexRefiner::refine_via_edges(), libMesh::SimplexRefiner::refine_via_edges(), libMesh::StaticCondensationDofMap::reinit(), libMesh::DistributedMesh::renumber_dof_objects(), libMesh::DistributedMesh::renumber_nodes_and_elements(), libMesh::DofMap::scatter_constraints(), libMesh::CheckpointIO::select_split_config(), libMesh::DistributedMesh::set_next_unique_id(), libMesh::DofMap::set_nonlocal_dof_objects(), libMesh::PetscDMWrapper::set_point_range_in_section(), libMesh::RBEIMEvaluation::side_gather_bfs(), MeshFunctionTest::test_bad_gradient_var_with_out_of_mesh_value(), MeshFunctionTest::test_bad_hessian_var_with_out_of_mesh_value(), ExodusTest< elem_type >::test_read_gold(), ExodusTest< elem_type >::test_write(), ExodusC0PolyhedronTest::test_write_and_read_hexagonal_prism(), ExodusC0PolygonTest::test_write_and_read_pentagon(), MeshInputTest::testAbaqusRead(), MeshInputTest::testBadGmsh(), BoundaryInfoTest::testBoundaryIDs(), MeshInputTest::testCopyElementSolutionImpl(), MeshInputTest::testCopyElementVectorImpl(), MeshInputTest::testCopyNodalSolutionImpl(), DefaultCouplingTest::testCoupling(), PointNeighborCouplingTest::testCoupling(), MeshInputTest::testDynaFileMappings(), MeshInputTest::testDynaNoSplines(), MeshInputTest::testDynaReadElem(), MeshInputTest::testDynaReadPatch(), MeshInputTest::testExodusFileMappings(), MeshInputTest::testExodusIGASidesets(), MeshInputTest::testExodusWriteElementDataFromDiscontinuousNodalData(), MeshInputTest::testGmshBCIDOverlap(), MeshInputTest::testGoodGmsh(), MeshInputTest::testGoodSTL(), MeshInputTest::testGoodSTLBinary(), BoundaryInfoTest::testInternalBoundary(), MeshInputTest::testLowOrderEdgeBlocks(), BoundaryMeshSubdomainTest::testPerBoundarySubdomain(), SystemsTest::testProjectMatrix3D(), BoundaryInfoTest::testShellFaceConstraints(), MeshInputTest::testSingleElementImpl(), BoundaryMeshSubdomainTest::testSingleSubdomain(), WriteVecAndScalar::testSolution(), CheckpointIOTest::testSplitter(), MeshInputTest::testTetgenIO(), MeshSmootherTest::testVariationalSmoother(), libMesh::MeshTools::total_weight(), libMesh::NetGenMeshInterface::triangulate(), libMesh::Parallel::Packing< Elem * >::unpack(), libMesh::Parallel::Packing< Node * >::unpack(), libMesh::DistributedMesh::update_parallel_id_counts(), libMesh::DTKAdapter::update_variable_values(), libMesh::MeshTools::volume(), libMesh::STLIO::write(), libMesh::XdrIO::write(), libMesh::NameBasedIO::write(), libMesh::CheckpointIO::write(), libMesh::EquationSystems::write(), libMesh::EquationSystems::write(), libMesh::GMVIO::write_discontinuous_gmv(), libMesh::ExodusII_IO::write_element_data(), libMesh::ExodusII_IO::write_element_data_from_discontinuous_nodal_data(), libMesh::ExodusII_IO_Helper::write_element_values(), libMesh::ExodusII_IO_Helper::write_element_values_element_major(), libMesh::ExodusII_IO_Helper::write_elements(), libMesh::ExodusII_IO_Helper::write_elemset_data(), libMesh::ExodusII_IO_Helper::write_elemsets(), libMesh::ExodusII_IO::write_global_data(), libMesh::ExodusII_IO_Helper::write_global_values(), libMesh::System::write_header(), libMesh::ExodusII_IO::write_information_records(), libMesh::ExodusII_IO_Helper::write_information_records(), libMesh::ExodusII_IO_Helper::write_nodal_coordinates(), libMesh::ExodusII_IO::write_nodal_data(), libMesh::VTKIO::write_nodal_data(), libMesh::UCDIO::write_nodal_data(), libMesh::ExodusII_IO::write_nodal_data_common(), libMesh::ExodusII_IO::write_nodal_data_discontinuous(), libMesh::ExodusII_IO_Helper::write_nodal_values(), libMesh::ExodusII_IO_Helper::write_nodeset_data(), libMesh::ExodusII_IO_Helper::write_nodesets(), libMesh::Nemesis_IO_Helper::write_nodesets(), libMesh::RBEIMEvaluation::write_out_interior_basis_functions(), libMesh::RBEIMEvaluation::write_out_node_basis_functions(), libMesh::RBEIMEvaluation::write_out_side_basis_functions(), write_output_solvedata(), libMesh::System::write_parallel_data(), libMesh::RBConstruction::write_riesz_representors_to_files(), libMesh::System::write_SCALAR_dofs(), libMesh::XdrIO::write_serialized_bc_names(), libMesh::XdrIO::write_serialized_bcs_helper(), libMesh::System::write_serialized_blocked_dof_objects(), libMesh::XdrIO::write_serialized_connectivity(), libMesh::System::write_serialized_data(), libMesh::XdrIO::write_serialized_nodes(), libMesh::XdrIO::write_serialized_nodesets(), libMesh::XdrIO::write_serialized_subdomain_names(), libMesh::System::write_serialized_vector(), libMesh::System::write_serialized_vectors(), libMesh::ExodusII_IO_Helper::write_sideset_data(), libMesh::ExodusII_IO_Helper::write_sidesets(), libMesh::Nemesis_IO_Helper::write_sidesets(), libMesh::ExodusII_IO::write_timestep(), libMesh::ExodusII_IO_Helper::write_timestep(), and libMesh::ExodusII_IO::write_timestep_discontinuous().

◆ read()

template<typename T >
void libMesh::SparseMatrix< T >::read ( const std::string &  filename)
virtualinherited

Read the contents of the matrix from a file, with the file format inferred from the extension of filename.

Definition at line 686 of file sparse_matrix.C.

687{
688 {
689 std::ifstream in (filename.c_str());
690 libmesh_error_msg_if
691 (!in.good(), "ERROR: cannot read file:\n\t" <<
692 filename);
693 }
694
695 std::string_view basename = Utility::basename_of(filename);
696
697 const bool gzipped_file = Utility::ends_with(filename, ".gz");
698
699 if (gzipped_file)
700 basename.remove_suffix(3);
701
702 if (Utility::ends_with(basename, ".matlab") ||
703 Utility::ends_with(basename, ".m"))
704 this->read_matlab(filename);
705 else if (Utility::ends_with(basename, ".petsc64"))
706 {
707#ifndef LIBMESH_HAVE_PETSC
708 libmesh_error_msg("Cannot load PETSc matrix file " <<
709 filename << " without PETSc-enabled libMesh.");
710#elif LIBMESH_DOF_ID_BYTES != 8
711 libmesh_error_msg("Cannot load 64-bit PETSc matrix file " <<
712 filename << " with non-64-bit libMesh.");
713#endif
714 if (gzipped_file)
715 libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
716 this->read_petsc_binary(filename);
717 }
718 else if (Utility::ends_with(basename, ".petsc32"))
719 {
720#ifndef LIBMESH_HAVE_PETSC
721 libmesh_error_msg("Cannot load PETSc matrix file " <<
722 filename << " without PETSc-enabled libMesh.");
723#elif LIBMESH_DOF_ID_BYTES != 4
724 libmesh_error_msg("Cannot load 32-bit PETSc matrix file " <<
725 filename << " with non-32-bit libMesh.");
726#endif
727 if (gzipped_file)
728 libmesh_not_implemented_msg("Gzipped PETSc matrices are not currently supported");
729 this->read_petsc_binary(filename);
730 }
731 else if (Utility::ends_with(basename, ".h5"))
732 {
733 if (gzipped_file)
734 libmesh_not_implemented_msg("Gzipped HDF5 matrices are not currently supported");
735 this->read_coreform_hdf5(filename);
736 }
737 else
738 libmesh_error_msg(" ERROR: Unrecognized matrix file extension on: "
739 << basename
740 << "\n I understand the following:\n\n"
741 << " *.h5 -- CoreForm HDF5 sparse matrix format\n"
742 << " *.matlab -- Matlab sparse matrix format\n"
743 << " *.matlab.gz -- Matlab sparse matrix format, gzipped\n"
744 << " *.m -- Matlab sparse matrix format\n"
745 << " *.m.gz -- Matlab sparse matrix format, gzipped\n"
746 << " *.petsc32 -- PETSc binary format, 32-bit\n"
747 << " *.petsc64 -- PETSc binary format, 64-bit\n"
748 );
749}
virtual void read_coreform_hdf5(const std::string &filename, const std::string &groupname="extraction")
Read the contents of the matrix from a file, with the HDF5 sparse matrix format used by CoreForm,...
virtual void read_matlab(const std::string &filename)
Read the contents of the matrix from the Matlab-script sparse matrix format used by PETSc.
virtual void read_petsc_binary(const std::string &filename)
Read the contents of the matrix from a file in PETSc's binary sparse matrix format.

References libMesh::Utility::basename_of(), and libMesh::Utility::ends_with().

◆ read_coreform_hdf5()

template<typename T >
void libMesh::SparseMatrix< T >::read_coreform_hdf5 ( const std::string &  filename,
const std::string &  groupname = "extraction" 
)
virtualinherited

Read the contents of the matrix from a file, with the HDF5 sparse matrix format used by CoreForm, expecing sparse matrix data in the group given by groupname.

This will be initialized with the sparsity from the file, linearly partitioned onto the number of processors available unless this matrix is pre-sized and pre-partitionsed.

Definition at line 822 of file sparse_matrix.C.

824{
825#ifndef LIBMESH_HAVE_HDF5
826 libmesh_ignore(filename, groupname);
827 libmesh_error_msg("ERROR: need HDF5 support to handle .h5 files!!!");
828#else
829 LOG_SCOPE("read_coreform_hdf5()", "SparseMatrix");
830
831 std::size_t num_rows = 0, num_cols = 0;
832
833 // These are only used on pid 0, but avoid "uninitialized" warnings
834 hid_t group = H5I_INVALID_HID;
835 hid_t file = H5I_INVALID_HID;
836
837 if (this->processor_id() == 0)
838 {
839 file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
840
841 if (file == H5I_INVALID_HID)
842 libmesh_file_error(filename);
843
844 group = H5Gopen(file, groupname.c_str(), H5P_DEFAULT);
845 check_open(filename, group, groupname);
846
847 auto read_size_attribute = [&filename, &group]
848 (const std::string & attribute_name)
849 {
850 unsigned long long returnval = 0;
851
852 const hid_t attr = H5Aopen(group, attribute_name.c_str(), H5P_DEFAULT);
853 check_open(filename, attr, attribute_name);
854
855 const hid_t attr_type = H5Aget_type(attr);
856 check_hdf5(filename, attr_type, attribute_name + " type");
857
858 // HDF5 will convert between the file's integer type and ours, but
859 // we do expect an integer type.
860 if (H5Tget_class(attr_type) != H5T_INTEGER)
861 libmesh_error_msg("Non-integer type for " + attribute_name + " in " + filename);
862
863 H5Tclose(attr_type);
864
865 // HDF5 is supposed to handle both upscaling and endianness
866 // conversions here
867 const herr_t errval = H5Aread(attr, H5T_NATIVE_ULLONG, &returnval);
868 check_hdf5(filename, errval, attribute_name + " read");
869
870 H5Aclose(attr);
871
872 return returnval;
873 };
874
875 num_cols = read_size_attribute("num_cols");
876 num_rows = read_size_attribute("num_rows");
877
878 this->comm().broadcast(num_cols);
879 this->comm().broadcast(num_rows);
880 }
881 else
882 {
883 this->comm().broadcast(num_cols);
884 this->comm().broadcast(num_rows);
885 }
886
887 numeric_index_type new_row_start, new_row_stop,
888 new_col_start, new_col_stop;
889
890 // If we need to reinit, we need to determine which rows+columns
891 // each processor is in charge of.
892 std::vector<numeric_index_type> new_row_starts, new_row_stops,
893 new_col_starts, new_col_stops;
894
895 if (this->initialized() &&
896 num_cols == this->n() &&
897 num_rows == this->m())
898 {
899 new_row_start = this->row_start(),
900 new_row_stop = this->row_stop();
901
902 new_col_start = this->col_start(),
903 new_col_stop = this->col_stop();
904 }
905 else
906 {
907 // Determine which rows/columns each processor will be in charge of
908 new_row_start = this->processor_id() * num_rows / this->n_processors(),
909 new_row_stop = (this->processor_id()+1) * num_rows / this->n_processors();
910
911 new_col_start = this->processor_id() * num_cols / this->n_processors(),
912 new_col_stop = (this->processor_id()+1) * num_cols / this->n_processors();
913 }
914
915 this->comm().gather(0, new_row_start, new_row_starts);
916 this->comm().gather(0, new_row_stop, new_row_stops);
917 this->comm().gather(0, new_col_start, new_col_starts);
918 this->comm().gather(0, new_col_stop, new_col_stops);
919
920 numeric_index_type on_diagonal_nonzeros = 0,
921 off_diagonal_nonzeros = 0;
922
923 std::vector<std::size_t> cols, row_offsets;
924 std::vector<double> vals;
925
926 if (this->processor_id() == 0)
927 {
928 auto read_vector = [&filename, &group]
929 (const std::string & dataname, auto hdf5_class,
930 auto hdf5_type, auto & datavec)
931 {
932 const hid_t data = H5Dopen1(group, dataname.c_str());
933 check_open(filename, data, dataname.c_str());
934
935 const hid_t data_type = H5Dget_type(data);
936 check_hdf5(filename, data_type, dataname + " type");
937
938 // HDF5 will convert between the file's integer type and ours, but
939 // we do expect an integer type.
940 if (H5Tget_class(data_type) != hdf5_class)
941 libmesh_error_msg("Unexpected type for " + dataname + " in " + filename);
942
943 H5Tclose(data_type);
944
945 const hid_t dataspace = H5Dget_space(data);
946 check_hdf5(filename, dataspace, dataname + " space");
947
948 int ndims = H5Sget_simple_extent_ndims(dataspace);
949 if (ndims != 1)
950 libmesh_error_msg("Non-vector space for " + dataname + " in " + filename);
951
952 hsize_t len, maxlen;
953 herr_t errval = H5Sget_simple_extent_dims(dataspace, &len, &maxlen);
954 check_hdf5(filename, errval, dataname + " dims");
955
956 datavec.resize(len);
957
958 errval = H5Dread(data, hdf5_type, H5S_ALL, H5S_ALL, H5P_DEFAULT, datavec.data());
959 check_hdf5(filename, errval, dataname + " read");
960
961 H5Dclose(data);
962 };
963
964 read_vector("cols", H5T_INTEGER, H5T_NATIVE_ULLONG, cols);
965 read_vector("row_offsets", H5T_INTEGER, H5T_NATIVE_ULLONG, row_offsets);
966 read_vector("vals", H5T_FLOAT, H5T_NATIVE_DOUBLE, vals);
967
968 if (cols.size() != vals.size())
969 libmesh_error_msg("Inconsistent cols/vals sizes in " + filename);
970
971 if (row_offsets.size() != num_rows + 1)
972 libmesh_error_msg("Inconsistent row_offsets size in " + filename);
973
974 // Data for the row we're working on
975 numeric_index_type current_row = 0;
976 processor_id_type current_proc = 0;
977 numeric_index_type current_on_diagonal_nonzeros = 0;
978 numeric_index_type current_off_diagonal_nonzeros = 0;
979 if (row_offsets[0] != 0)
980 libmesh_error_msg("Unexpected row_offsets[0] in " + filename);
981
982 for (auto i : index_range(cols))
983 {
984 while (row_offsets[current_row+1] <= i)
985 {
986 ++current_row;
987 if (row_offsets[current_row] < row_offsets[current_row-1])
988 libmesh_error_msg("Non-monotonic row_offsets in " + filename);
989 current_on_diagonal_nonzeros = 0;
990 current_off_diagonal_nonzeros = 0;
991 }
992
993 while (current_row >= new_row_stops[current_proc])
994 ++current_proc;
995
996 // 0-based indexing in file
997 if (cols[i] >= new_col_starts[current_proc] &&
998 cols[i] < new_col_stops[current_proc])
999 {
1000 ++current_on_diagonal_nonzeros;
1001 on_diagonal_nonzeros =
1002 std::max(on_diagonal_nonzeros,
1003 current_on_diagonal_nonzeros);
1004 }
1005 else
1006 {
1007 ++current_off_diagonal_nonzeros;
1008 off_diagonal_nonzeros =
1009 std::max(off_diagonal_nonzeros,
1010 current_off_diagonal_nonzeros);
1011 }
1012 }
1013 }
1014
1015 this->comm().broadcast(on_diagonal_nonzeros);
1016 this->comm().broadcast(off_diagonal_nonzeros);
1017
1018 this->init(num_rows, num_cols,
1019 new_row_stop-new_row_start,
1020 new_col_stop-new_col_start,
1021 on_diagonal_nonzeros,
1022 off_diagonal_nonzeros);
1023
1024 // Set the matrix values last.
1025 if (this->processor_id() == 0)
1026 {
1027 numeric_index_type current_row = 0;
1028 for (auto i : index_range(cols))
1029 {
1030 while (row_offsets[current_row+1] <= i)
1031 {
1032 ++current_row;
1033 libmesh_assert_greater_equal (row_offsets[current_row],
1034 row_offsets[current_row-1]);
1035 }
1036 this->set(current_row, cols[i], vals[i]);
1037 }
1038
1039 H5Gclose(group);
1040 H5Fclose(file);
1041 }
1042
1043 this->close();
1044#endif // LIBMESH_HAVE_HDF5
1045}
virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1)=0
Initialize SparseMatrix with the specified sizes.
virtual void close()=0
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
virtual numeric_index_type col_stop() const =0
virtual numeric_index_type col_start() const =0
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value)=0
Set the element (i,j) to value.
MPI_Datatype data_type

References libMesh::index_range(), libMesh::initialized(), and libMesh::libmesh_ignore().

◆ read_matlab()

template<typename T >
void libMesh::SparseMatrix< T >::read_matlab ( const std::string &  filename)
virtualinherited

Read the contents of the matrix from the Matlab-script sparse matrix format used by PETSc.

If the size and sparsity of the matrix in filename appears consistent with the existing sparsity of this then the existing parallel decomposition and sparsity will be retained. If not, then this will be initialized with the sparsity from the file, linearly partitioned onto the number of processors available.

Definition at line 1050 of file sparse_matrix.C.

1051{
1052 LOG_SCOPE("read_matlab()", "SparseMatrix");
1053
1054#ifndef LIBMESH_HAVE_CXX11_REGEX
1055 libmesh_not_implemented(); // What is your compiler?!? Email us!
1056 libmesh_ignore(filename);
1057#else
1058 parallel_object_only();
1059
1060 const bool gzipped_file = Utility::ends_with(filename, ".gz");
1061
1062 // The sizes we get from the file
1063 std::size_t m = 0,
1064 n = 0;
1065
1066 // If we don't already have this size, we'll need to reinit, and
1067 // determine which rows+columns each processor is in charge of.
1068 std::vector<numeric_index_type> new_row_starts, new_row_stops,
1069 new_col_starts, new_col_stops;
1070
1071 numeric_index_type new_row_start, new_row_stop,
1072 new_col_start, new_col_stop;
1073
1074 // We'll read through the file three times: once to get a reliable
1075 // value for the matrix size (so we can divvy it up among
1076 // processors), then again to get the sparsity to send to each
1077 // processor, then a final time to get the entries to send to each
1078 // processor.
1079 //
1080 // We'll use an istream here; it might be an ifstream if we're
1081 // opening a raw ASCII file or a gzstream if we're opening a
1082 // compressed one.
1083 std::unique_ptr<std::istream> file;
1084
1085 // We'll need a temporary structure to cache matrix entries, because
1086 // we need to read through the whole file before we know the size
1087 // and sparsity structure with which we can init().
1088 //
1089 // Reading through the file three times via `seekg` doesn't work
1090 // with our gzstream wrapper, and seems to take three times as long
1091 // even with a plain ifstream. What happened to disk caching!?
1092 std::vector<std::tuple<numeric_index_type, numeric_index_type, T>> entries;
1093
1094 // First read through the file, saving size and entry data
1095 {
1096 // We'll read the matrix on processor 0 rather than try to juggle
1097 // parallel I/O.
1098 if (this->processor_id() == 0)
1099 {
1100 // We'll be using regular expressions to make ourselves slightly
1101 // more robust to formatting.
1102 const std::regex start_regex // assignment like "zzz = ["
1103 ("\\s*\\w+\\s*=\\s*\\[");
1104 const std::regex end_regex // end of assignment
1105 ("^[^%]*\\]");
1106
1107 if (gzipped_file)
1108 {
1109#ifdef LIBMESH_HAVE_GZSTREAM
1110 auto inf = std::make_unique<igzstream>();
1111 libmesh_assert(inf);
1112 inf->open(filename.c_str(), std::ios::in);
1113 file = std::move(inf);
1114#else
1115 libmesh_error_msg("ERROR: need gzstream to handle .gz files!!!");
1116#endif
1117 }
1118 else
1119 {
1120 auto inf = std::make_unique<std::ifstream>();
1121 libmesh_assert(inf);
1122
1123 std::string new_name = Utility::unzip_file(filename);
1124
1125 inf->open(new_name.c_str(), std::ios::in);
1126 file = std::move(inf);
1127 }
1128
1129 // If we have a matrix with all-zero trailing rows, the only
1130 // way to get the size is if it ended up in a comment
1131 const std::regex size_regex // comment like "% size = 8 8"
1132 ("%\\s*[Ss][Ii][Zz][Ee]\\s*=\\s*(\\d+)\\s+(\\d+)");
1133 const std::string whitespace = " \t";
1134
1135 bool have_started = false;
1136 bool have_ended = false;
1137 std::size_t largest_i_seen = 0, largest_j_seen = 0;
1138
1139 // Data for the row we're working on
1140 // Use 1-based indexing for current_row, as in the file
1141 std::size_t current_row = 1;
1142
1143 for (std::string line; std::getline(*file, line);)
1144 {
1145 std::smatch sm;
1146
1147 // First, try to match an entry. This is the most common
1148 // case so we won't rely on slow std::regex for it.
1149 // stringstream is at least an improvement over that.
1150
1151 // Look for row/col/val like "1 1 -2.0e-4"
1152
1153 std::istringstream l(line);
1154
1155 std::size_t i, j;
1156 T value;
1157
1158 l >> i >> j >> value;
1159
1160 if (!l.fail())
1161 {
1162 libmesh_error_msg_if
1163 (!have_started, "Confused by premature entries in matrix file " << filename);
1164
1165 entries.emplace_back(cast_int<numeric_index_type>(i),
1166 cast_int<numeric_index_type>(j),
1167 value);
1168
1169 libmesh_error_msg_if
1170 (!i || !j, "Expected 1-based indexing in matrix file "
1171 << filename);
1172
1173 current_row = std::max(current_row, i);
1174
1175 libmesh_error_msg_if
1176 (i < current_row,
1177 "Can't handle out-of-order entries in matrix file "
1178 << filename);
1179
1180 largest_i_seen = std::max(i, largest_i_seen);
1181 largest_j_seen = std::max(j, largest_j_seen);
1182 }
1183
1184 else if (std::regex_search(line, sm, size_regex))
1185 {
1186 const std::string msize = sm[1];
1187 const std::string nsize = sm[2];
1188 m = std::stoull(msize);
1189 n = std::stoull(nsize);
1190 }
1191
1192 else if (std::regex_search(line, start_regex))
1193 have_started = true;
1194
1195 else if (std::regex_search(line, end_regex))
1196 {
1197 have_ended = true;
1198 break;
1199 }
1200 }
1201
1202 libmesh_error_msg_if
1203 (!have_started, "Confused by missing assignment beginning in matrix file " << filename);
1204
1205 libmesh_error_msg_if
1206 (!have_ended, "Confused by missing assignment ending in matrix file " << filename);
1207
1208 libmesh_error_msg_if
1209 (m > largest_i_seen, "Confused by missing final row(s) in matrix file " << filename);
1210
1211 libmesh_error_msg_if
1212 (m > 0 && m < largest_i_seen, "Confused by extra final row(s) in matrix file " << filename);
1213
1214 if (!m)
1215 m = largest_i_seen;
1216
1217 libmesh_error_msg_if
1218 (n > largest_j_seen, "Confused by missing final column(s) in matrix file " << filename);
1219
1220 libmesh_error_msg_if
1221 (n > 0 && n < largest_j_seen, "Confused by extra final column(s) in matrix file " << filename);
1222
1223 if (!n)
1224 n = largest_j_seen;
1225
1226 this->comm().broadcast(m);
1227 this->comm().broadcast(n);
1228 }
1229 else
1230 {
1231 this->comm().broadcast(m);
1232 this->comm().broadcast(n);
1233 }
1234
1235 if (this->initialized() &&
1236 m == this->m() &&
1237 n == this->n())
1238 {
1239 new_row_start = this->row_start(),
1240 new_row_stop = this->row_stop();
1241
1242 new_col_start = this->col_start(),
1243 new_col_stop = this->col_stop();
1244 }
1245 else
1246 {
1247 // Determine which rows/columns each processor will be in charge of
1248 new_row_start = this->processor_id() * m / this->n_processors(),
1249 new_row_stop = (this->processor_id()+1) * m / this->n_processors();
1250
1251 new_col_start = this->processor_id() * n / this->n_processors(),
1252 new_col_stop = (this->processor_id()+1) * n / this->n_processors();
1253 }
1254
1255 this->comm().gather(0, new_row_start, new_row_starts);
1256 this->comm().gather(0, new_row_stop, new_row_stops);
1257 this->comm().gather(0, new_col_start, new_col_starts);
1258 this->comm().gather(0, new_col_stop, new_col_stops);
1259
1260 } // Done reading entry data and broadcasting matrix size
1261
1262 // Calculate the matrix sparsity and initialize it second
1263 {
1264 // Deduce the sparsity pattern, or at least the maximum number of
1265 // on- and off- diagonal non-zeros per row.
1266 numeric_index_type on_diagonal_nonzeros =0,
1267 off_diagonal_nonzeros =0;
1268
1269 if (this->processor_id() == 0)
1270 {
1271 // Data for the row we're working on
1272 // Use 1-based indexing for current_row, as in the file
1273 numeric_index_type current_row = 1;
1274 processor_id_type current_proc = 0;
1275 numeric_index_type current_on_diagonal_nonzeros = 0;
1276 numeric_index_type current_off_diagonal_nonzeros = 0;
1277
1278 for (auto [i, j, value] : entries)
1279 {
1280 if (i > current_row)
1281 {
1282 current_row = i;
1283 // +1 for 1-based indexing in file
1284 while (current_row >= (new_row_stops[current_proc]+1))
1285 ++current_proc;
1286 current_on_diagonal_nonzeros = 0;
1287 current_off_diagonal_nonzeros = 0;
1288 }
1289
1290 // +1 for 1-based indexing in file
1291 if (j >= (new_col_starts[current_proc]+1) &&
1292 j < (new_col_stops[current_proc]+1))
1293 {
1294 ++current_on_diagonal_nonzeros;
1295 on_diagonal_nonzeros =
1296 std::max(on_diagonal_nonzeros,
1297 current_on_diagonal_nonzeros);
1298 }
1299 else
1300 {
1301 ++current_off_diagonal_nonzeros;
1302 off_diagonal_nonzeros =
1303 std::max(off_diagonal_nonzeros,
1304 current_off_diagonal_nonzeros);
1305 }
1306 }
1307 }
1308
1309 this->comm().broadcast(on_diagonal_nonzeros);
1310 this->comm().broadcast(off_diagonal_nonzeros);
1311
1312 this->init(m, n,
1313 new_row_stop-new_row_start,
1314 new_col_stop-new_col_start,
1315 on_diagonal_nonzeros,
1316 off_diagonal_nonzeros);
1317 }
1318
1319 // Set the matrix values last.
1320 // Convert from 1-based to 0-based indexing
1321 if (this->processor_id() == 0)
1322 for (auto [i, j, value] : entries)
1323 this->set(i-1, j-1, value);
1324
1325 this->close();
1326#endif
1327}
std::string unzip_file(std::string_view name)
Create an unzipped copy of a bz2 or xz file, returning the name of the now-unzipped file that can be ...
Definition utility.C:164

References libMesh::Utility::ends_with(), libMesh::initialized(), libMesh::libmesh_assert(), libMesh::libmesh_ignore(), libMesh::Utility::unzip_file(), and value.

Referenced by ConstraintOperatorTest::test1DCoarseningNewNodes(), and ConstraintOperatorTest::testCoreform().

◆ read_petsc_binary()

template<typename T >
void libMesh::PetscMatrix< T >::read_petsc_binary ( const std::string &  filename)
overridevirtual

Read the contents of the matrix from a file in PETSc's binary sparse matrix format.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 715 of file petsc_matrix.C.

716{
717 LOG_SCOPE("read_petsc_binary()", "PetscMatrix");
718
719 this->_petsc_viewer(filename, PETSCVIEWERBINARY, FILE_MODE_READ);
720}

◆ read_petsc_hdf5()

template<typename T >
void libMesh::PetscMatrix< T >::read_petsc_hdf5 ( const std::string &  filename)
overridevirtual

Read the contents of the matrix from a file in PETSc's HDF5 sparse matrix format.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 725 of file petsc_matrix.C.

726{
727 LOG_SCOPE("read_petsc_hdf5()", "PetscMatrix");
728
729 this->_petsc_viewer(filename, PETSCVIEWERHDF5, FILE_MODE_READ);
730}

◆ reinit_submatrix()

template<typename T >
virtual void libMesh::SparseMatrix< T >::reinit_submatrix ( SparseMatrix< T > &  submatrix,
const std::vector< numeric_index_type > &  rows,
const std::vector< numeric_index_type > &  cols 
) const
inlinevirtualinherited

This function is similar to the one above, but it allows you to reuse the existing sparsity pattern of "submatrix" instead of reallocating it again.

This should hopefully be more efficient if you are frequently extracting submatrices of the same size.

Definition at line 565 of file sparse_matrix.h.

568 {
569 this->_get_submatrix(submatrix,
570 rows,
571 cols,
572 true); // true means REUSE submatrix
573 }

References libMesh::SparseMatrix< T >::_get_submatrix().

◆ require_sparsity_pattern()

template<typename T >
virtual bool libMesh::SparseMatrix< T >::require_sparsity_pattern ( ) const
inlinevirtualinherited
Returns
Whether this matrix needs the sparsity pattern computed by the DofMap

Reimplemented in libMesh::PetscMatrixShellMatrix< T >, libMesh::PetscMatrixShellMatrix< Number >, and libMesh::StaticCondensation.

Definition at line 168 of file sparse_matrix.h.

168{ return !this->use_hash_table(); }
bool use_hash_table() const

References libMesh::SparseMatrix< T >::use_hash_table().

◆ reset_preallocation()

template<typename T >
void libMesh::PetscMatrix< T >::reset_preallocation ( )

Reset matrix to use the original nonzero pattern provided by users.

Definition at line 385 of file petsc_matrix.C.

386{
387 semiparallel_only();
388
389#if !PETSC_VERSION_LESS_THAN(3,9,0)
390 libmesh_assert (this->initialized());
391
392 LibmeshPetscCall(MatResetPreallocation(this->_mat));
393#else
394 libmesh_warning("Your version of PETSc doesn't support resetting of "
395 "preallocation, so we will use your most recent sparsity "
396 "pattern. This may result in a degradation of performance\n");
397#endif
398}

References libMesh::initialized(), and libMesh::libmesh_assert().

◆ restore_original_nonzero_pattern()

template<typename T >
void libMesh::PetscMatrix< T >::restore_original_nonzero_pattern ( )
overridevirtual

Reset the memory storage of the matrix.

Unlike 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 use_hash_table()) will reset (clear) the hash table used for assembly. In the words of the MatResetPreallocation documentation in PETSc, 'current values in the matrix are lost in this call', so a user can expect to have back their original sparsity pattern in a zeroed state

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1267 of file petsc_matrix.C.

1268{
1269 semiparallel_only();
1270
1271 if (this->_use_hash_table)
1272#if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
1273 // This performs MatReset plus re-establishes the hash table
1274 LibmeshPetscCall(MatResetHash(this->_mat));
1275#else
1276 libmesh_error_msg("Resetting hash tables not supported until PETSc version 3.23");
1277#endif
1278 else
1279 this->reset_preallocation();
1280}
void reset_preallocation()
Reset matrix to use the original nonzero pattern provided by users.
bool _use_hash_table
Flag indicating whether the matrix is assembled using a hash table.

◆ row_start()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::row_start ( ) const
overridevirtualinherited
Returns
The index of the first matrix row stored on this processor.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 178 of file petsc_matrix_base.C.

179{
180 libmesh_assert (this->initialized());
181
182 PetscInt start=0, stop=0;
183
184 LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
185
186 return static_cast<numeric_index_type>(start);
187}

References libMesh::initialized(), and libMesh::libmesh_assert().

Referenced by libMesh::PetscPreconditioner< T >::set_hypre_ads_data(), and libMesh::PetscPreconditioner< T >::set_hypre_ams_data().

◆ row_stop()

template<typename T >
numeric_index_type libMesh::PetscMatrixBase< T >::row_stop ( ) const
overridevirtualinherited
Returns
The index of the last matrix row (+1) stored on this processor.

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 190 of file petsc_matrix_base.C.

191{
192 libmesh_assert (this->initialized());
193
194 PetscInt start=0, stop=0;
195
196 LibmeshPetscCall(MatGetOwnershipRange(this->_mat, &start, &stop));
197
198 return static_cast<numeric_index_type>(stop);
199}

References libMesh::initialized(), and libMesh::libmesh_assert().

Referenced by libMesh::PetscPreconditioner< T >::set_hypre_ads_data(), and libMesh::PetscPreconditioner< T >::set_hypre_ams_data().

◆ scale()

template<typename T >
void libMesh::PetscMatrix< T >::scale ( const T  scale)
overridevirtual

Scales all elements of this matrix by scale.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1234 of file petsc_matrix.C.

1235{
1236 libmesh_assert(this->closed());
1237
1238 LibmeshPetscCall(MatScale(this->_mat, PS(scale)));
1239}
virtual void scale(const T scale) override
Scales all elements of this matrix by scale.

References libMesh::closed(), libMesh::libmesh_assert(), and libMesh::PS().

◆ set()

template<typename T >
void libMesh::PetscMatrix< T >::set ( const numeric_index_type  i,
const numeric_index_type  j,
const T  value 
)
overridevirtual

Set the element (i,j) to value.

Throws an error if the entry does not exist. Zero values can be "stored" in non-existent fields.

Implements libMesh::SparseMatrix< T >.

Definition at line 972 of file petsc_matrix.C.

975{
976 libmesh_assert (this->initialized());
977
978 PetscInt i_val=i, j_val=j;
979
980 PetscScalar petsc_value = static_cast<PetscScalar>(value);
981 std::scoped_lock lock(this->_petsc_matrix_mutex);
982 LibmeshPetscCall(MatSetValues(this->_mat, 1, &i_val, 1, &j_val,
983 &petsc_value, INSERT_VALUES));
984}

References libMesh::initialized(), libMesh::libmesh_assert(), and value.

Referenced by libMesh::PetscPreconditioner< T >::set_hypre_ads_data(), and libMesh::PetscPreconditioner< T >::set_hypre_ams_data().

◆ set_context()

template<typename T >
void libMesh::PetscMatrixBase< T >::set_context ( )
inherited

Set the context (ourself) for _mat.

Definition at line 105 of file petsc_matrix_base.C.

106{
107 libmesh_assert(this->_mat);
108 PetscContainer container;
109 LibmeshPetscCall(PetscContainerCreate(this->comm().get(), &container));
110 LibmeshPetscCall(PetscContainerSetPointer(container, this));
111 LibmeshPetscCall(PetscObjectCompose((PetscObject)(Mat)this->_mat, "PetscMatrixCtx", (PetscObject)container));
112 LibmeshPetscCall(PetscContainerDestroy(&container));
113}

References libMesh::libmesh_assert().

Referenced by libMesh::PetscMatrixBase< T >::PetscMatrixBase().

◆ set_destroy_mat_on_exit()

template<typename T >
void libMesh::PetscMatrixBase< T >::set_destroy_mat_on_exit ( bool  destroy = true)
inherited

If set to false, we don't delete the Mat on destruction and allow instead for PETSc to manage it.

Definition at line 91 of file petsc_matrix_base.C.

92{
94}
void destroy(triangulateio &t, IO_Type)
Frees any memory which has been dynamically allocated by Triangle.

◆ solver_package()

template<typename T >
virtual SolverPackage libMesh::PetscMatrixBase< T >::solver_package ( )
inlineoverridevirtualinherited

Implements libMesh::SparseMatrix< T >.

Reimplemented in libMesh::StaticCondensation, and libMesh::StaticCondensation.

Definition at line 107 of file petsc_matrix_base.h.

108 {
109 return PETSC_SOLVERS;
110 }

References libMesh::PETSC_SOLVERS.

◆ supports_hash_table()

template<typename T >
bool libMesh::PetscMatrix< T >::supports_hash_table ( ) const
overridevirtual
Returns
Whether the matrix supports hash table assembly

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 1242 of file petsc_matrix.C.

1243{
1244#if PETSC_RELEASE_LESS_THAN(3,19,0)
1245 return false;
1246#else
1247 return true;
1248#endif
1249}

◆ swap()

template<typename T >
void libMesh::PetscMatrixBase< T >::swap ( PetscMatrixBase< T > &  m_in)
inherited

Swaps the internal data pointers of two PetscMatrices, no actual values are swapped.

Definition at line 98 of file petsc_matrix_base.C.

99{
100 std::swap(_mat, m_in._mat);
101 std::swap(_destroy_mat_on_exit, m_in._destroy_mat_on_exit);
102}

References libMesh::PetscMatrixBase< T >::_destroy_mat_on_exit, and libMesh::PetscMatrixBase< T >::_mat.

Referenced by DMlibMeshJacobian().

◆ update_preallocation_and_zero()

template<typename T >
void libMesh::PetscMatrix< T >::update_preallocation_and_zero ( )

Update the sparsity pattern based on dof_map, and set the matrix to zero.

This is useful in cases where the sparsity pattern changes during a computation.

Definition at line 379 of file petsc_matrix.C.

380{
381 libmesh_not_implemented();
382}

◆ update_sparsity_pattern()

template<typename T >
virtual void libMesh::SparseMatrix< T >::update_sparsity_pattern ( const SparsityPattern::Graph )
inlinevirtualinherited

Updates the matrix sparsity pattern.

When your SparseMatrix<T> implementation does not need this data, simply do not override this method.

Reimplemented in libMesh::LaspackMatrix< T >, and libMesh::EpetraMatrix< T >.

Definition at line 175 of file sparse_matrix.h.

175{}

Referenced by libMesh::DofMap::update_sparsity_pattern().

◆ use_hash_table() [1/2]

template<typename T >
bool libMesh::SparseMatrix< T >::use_hash_table ( ) const
inlineinherited
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

Definition at line 634 of file sparse_matrix.h.

634{ return _use_hash_table; }

References libMesh::SparseMatrix< T >::_use_hash_table.

Referenced by libMesh::SparseMatrix< T >::require_sparsity_pattern().

◆ use_hash_table() [2/2]

template<typename T >
void libMesh::SparseMatrix< T >::use_hash_table ( bool  use_hash)
inherited

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

Definition at line 692 of file sparse_matrix.h.

693{
694 libmesh_error_msg_if(use_hash && !this->supports_hash_table(),
695 "This matrix class does not support hash table assembly");
696 this->_use_hash_table = use_hash;
697}
virtual bool supports_hash_table() const

Referenced by PetscMatrixTest::testPetscCopyFromHash().

◆ vector_mult()

template<typename T >
void libMesh::SparseMatrix< T >::vector_mult ( NumericVector< T > &  dest,
const NumericVector< T > &  arg 
) const
inherited

Multiplies the matrix by the NumericVector arg and stores the result in NumericVector dest.

Definition at line 237 of file sparse_matrix.C.

239{
240 dest.zero();
241 this->vector_mult_add(dest,arg);
242}
void vector_mult_add(NumericVector< T > &dest, const NumericVector< T > &arg) const
Multiplies the matrix by the NumericVector arg and adds the result to the NumericVector dest.

References libMesh::NumericVector< T >::zero().

Referenced by libMesh::TransientRBConstruction::add_IC_to_RB_space(), libMesh::RBSCMConstruction::Aq_inner_product(), libMesh::AdvectionSystem::assemble_claw_rhs(), libMesh::RBSCMConstruction::B_inner_product(), libMesh::RBConstruction::compute_Fq_representor_innerprods(), libMesh::RBConstruction::compute_output_dual_innerprods(), libMesh::RBConstruction::compute_residual_dual_norm_slow(), libMesh::RBConstruction::enrich_RB_space(), libMesh::TransientRBConstruction::enrich_RB_space(), AssembleOptimization::gradient(), libMesh::TransientRBConstruction::mass_matrix_scaled_matvec(), AssembleOptimization::objective(), libMesh::RBConstruction::print_basis_function_orthogonality(), libMesh::ImplicitSystem::qoi_parameter_hessian(), libMesh::ImplicitSystem::qoi_parameter_hessian_vector_product(), libMesh::TransientRBConstruction::set_error_temporal_data(), libMesh::RBConstruction::train_reduced_basis_with_POD(), libMesh::TransientRBConstruction::truth_assembly(), libMesh::RBConstruction::truth_solve(), libMesh::RBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_RB_system_matrices(), libMesh::TransientRBConstruction::update_residual_terms(), and libMesh::RBConstruction::update_residual_terms().

◆ vector_mult_add()

template<typename T >
void libMesh::SparseMatrix< T >::vector_mult_add ( NumericVector< T > &  dest,
const NumericVector< T > &  arg 
) const
inherited

Multiplies the matrix by the NumericVector arg and adds the result to the NumericVector dest.

Definition at line 247 of file sparse_matrix.C.

249{
250 /* This functionality is actually implemented in the \p
251 NumericVector class. */
252 dest.add_vector(arg,*this);
253}

References libMesh::NumericVector< T >::add_vector().

Referenced by libMesh::ImplicitSystem::weighted_sensitivity_adjoint_solve().

◆ zero()

template<typename T >
void libMesh::PetscMatrix< T >::zero ( )
overridevirtual

Set all entries to 0.

Implements libMesh::SparseMatrix< T >.

Definition at line 401 of file petsc_matrix.C.

402{
403 libmesh_assert (this->initialized());
404
405 semiparallel_only();
406
407 PetscInt m_l, n_l;
408
409 LibmeshPetscCall(MatGetLocalSize(this->_mat,&m_l,&n_l));
410
411 if (n_l)
412 LibmeshPetscCall(MatZeroEntries(this->_mat));
413}

References libMesh::initialized(), and libMesh::libmesh_assert().

◆ zero_clone()

template<typename T >
std::unique_ptr< SparseMatrix< T > > libMesh::PetscMatrix< T >::zero_clone ( ) const
overridevirtual
Returns
A smart pointer to a copy of this matrix with the same type, size, and partitioning, but with all zero entries.
Note
This must be overridden in the derived classes.

Implements libMesh::SparseMatrix< T >.

Definition at line 435 of file petsc_matrix.C.

436{
437 libmesh_error_msg_if(!this->closed(), "Matrix must be closed before it can be cloned!");
438
439 // Copy the nonzero pattern only
440 Mat copy;
441 LibmeshPetscCall(MatDuplicate(this->_mat, MAT_DO_NOT_COPY_VALUES, &copy));
442
443 // Call wrapping PetscMatrix constructor, have it take over
444 // ownership.
445 auto ret = std::make_unique<PetscMatrix<T>>(copy, this->comm());
446 ret->set_destroy_mat_on_exit(true);
447
448 return ret;
449}

References libMesh::closed().

◆ zero_rows()

template<typename T >
void libMesh::PetscMatrix< T >::zero_rows ( std::vector< numeric_index_type > &  rows,
diag_value = 0.0 
)
overridevirtual

Sets all row entries to 0 then puts diag_value in the diagonal entry.

Reimplemented from libMesh::SparseMatrix< T >.

Definition at line 416 of file petsc_matrix.C.

417{
418 libmesh_assert (this->initialized());
419
420 semiparallel_only();
421
422 // As of petsc-dev at the time of 3.1.0, MatZeroRows now takes two additional
423 // optional arguments. The optional arguments (x,b) can be used to specify the
424 // solutions for the zeroed rows (x) and right hand side (b) to update.
425 // Could be useful for setting boundary conditions...
426 if (!rows.empty())
427 LibmeshPetscCall(MatZeroRows(this->_mat, cast_int<PetscInt>(rows.size()),
428 numeric_petsc_cast(rows.data()), PS(diag_value),
429 NULL, NULL));
430 else
431 LibmeshPetscCall(MatZeroRows(this->_mat, 0, NULL, PS(diag_value), NULL, NULL));
432}

References libMesh::initialized(), libMesh::libmesh_assert(), libMesh::numeric_petsc_cast(), and libMesh::PS().

Friends And Related Symbol Documentation

◆ ::PetscMatrixTest

template<typename T >
friend class ::PetscMatrixTest
friend

Definition at line 366 of file petsc_matrix.h.

Member Data Documentation

◆ _communicator

const Parallel::Communicator& libMesh::ParallelObject::_communicator
protectedinherited

◆ _counts

ReferenceCounter::Counts libMesh::ReferenceCounter::_counts
staticprotectedinherited

Actually holds the data.

Definition at line 124 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::get_info().

◆ _destroy_mat_on_exit

template<typename T >
bool libMesh::PetscMatrixBase< T >::_destroy_mat_on_exit
protectedinherited

This boolean value should only be set to false for the constructor which takes a PETSc Mat object.

Definition at line 186 of file petsc_matrix_base.h.

Referenced by libMesh::PetscMatrixBase< T >::swap().

◆ _dof_map

template<typename T >
DofMap const* libMesh::SparseMatrix< T >::_dof_map
protectedinherited

The DofMap object associated with this object.

May be queried for degree-of-freedom counts on processors.

Definition at line 666 of file sparse_matrix.h.

◆ _enable_print_counter

bool libMesh::ReferenceCounter::_enable_print_counter = true
staticprotectedinherited

Flag to control whether reference count information is printed when print_info is called.

Definition at line 143 of file reference_counter.h.

Referenced by libMesh::ReferenceCounter::disable_print_counter_info(), libMesh::ReferenceCounter::enable_print_counter_info(), and libMesh::ReferenceCounter::print_info().

◆ _is_initialized

template<typename T >
bool libMesh::SparseMatrix< T >::_is_initialized
protectedinherited

◆ _mat

template<typename T >
Mat libMesh::PetscMatrixBase< T >::_mat
protectedinherited

◆ _mat_type

template<typename T >
PetscMatrixType libMesh::PetscMatrix< T >::_mat_type
protected

Definition at line 349 of file petsc_matrix.h.

Referenced by libMesh::PetscMatrix< T >::PetscMatrix().

◆ _mutex

Threads::spin_mutex libMesh::ReferenceCounter::_mutex
staticprotectedinherited

Mutual exclusion object to enable thread-safe reference counting.

Definition at line 137 of file reference_counter.h.

◆ _n_objects

Threads::atomic< unsigned int > libMesh::ReferenceCounter::_n_objects
staticprotectedinherited

◆ _petsc_matrix_mutex [1/2]

template<typename T >
std::mutex libMesh::PetscMatrix< T >::_petsc_matrix_mutex
mutableprivate

Definition at line 353 of file petsc_matrix.h.

◆ _petsc_matrix_mutex [2/2]

template<typename T >
Threads::spin_mutex libMesh::PetscMatrix< T >::_petsc_matrix_mutex
mutableprivate

Definition at line 355 of file petsc_matrix.h.

◆ _sp

template<typename T >
SparsityPattern::Build const* libMesh::SparseMatrix< T >::_sp
protectedinherited

The sparsity pattern associated with this object.

Should be queried for entry counts (or with need_full_sparsity_pattern, patterns) when needed.

Definition at line 673 of file sparse_matrix.h.

◆ _use_hash_table

template<typename T >
bool libMesh::SparseMatrix< T >::_use_hash_table
protectedinherited

Flag indicating whether the matrix is assembled using a hash table.

Definition at line 683 of file sparse_matrix.h.

Referenced by libMesh::SparseMatrix< T >::use_hash_table().


The documentation for this class was generated from the following files: