12#ifndef MFEM_SPARSEMAT_HPP
13#define MFEM_SPARSEMAT_HPP
25#if defined(MFEM_USE_HIP)
26#if (HIP_VERSION_MAJOR * 100 + HIP_VERSION_MINOR) < 502
29#include <hipsparse/hipsparse.h>
38#if defined(__alignas_is_defined)
82#ifdef MFEM_USE_MEMALLOC
98#ifdef MFEM_USE_CUDA_OR_HIP
105#if defined(MFEM_USE_CUDA)
110#if CUDA_VERSION >= 10010
155 bool owna,
bool issorted);
224 inline const int *
GetI()
const {
return I; }
229 inline const int *
GetJ()
const {
return J; }
239 const int *
ReadI(
bool on_dev =
true)
const
255 const int *
ReadJ(
bool on_dev =
true)
const
285 int RowSize(
const int i)
const;
325 const real_t &
Elem(
int i,
int j)
const override;
353 const real_t a = 1.0)
const override;
366 const real_t a = 1.0)
const override;
527 real_t sc,
bool use_abs_diag =
false)
const;
531 real_t sc = 1.0,
bool use_abs_diag =
false)
const;
549 void Finalize(
int skip_zeros,
bool fix_empty_rows);
574 inline void SetColPtr(
const int row)
const;
599 void Set(
const int i,
const int j,
const real_t val);
600 void Add(
const int i,
const int j,
const real_t val);
740SparseMatrix *
Transpose(
const SparseMatrix &A);
751SparseMatrix *
Mult(
const SparseMatrix &A,
const SparseMatrix &B,
752 SparseMatrix *OAB = NULL);
755SparseMatrix *
TransposeMult(
const SparseMatrix &A,
const SparseMatrix &B);
759 const AbstractSparseMatrix &B);
762DenseMatrix *
Mult(
const SparseMatrix &A, DenseMatrix &B);
765DenseMatrix *
RAP(
const SparseMatrix &A, DenseMatrix &P);
768DenseMatrix *
RAP(DenseMatrix &A,
const SparseMatrix &P);
772SparseMatrix *
RAP(
const SparseMatrix &A,
const SparseMatrix &R,
773 SparseMatrix *ORAP = NULL);
776SparseMatrix *
RAP(
const SparseMatrix &Rt,
const SparseMatrix &A,
777 const SparseMatrix &P);
780SparseMatrix *
Mult_AtDA(
const SparseMatrix &A,
const Vector &D,
781 SparseMatrix *OAtDA = NULL);
785SparseMatrix *
Add(
const SparseMatrix & A,
const SparseMatrix & B);
788 const SparseMatrix & B);
790SparseMatrix *
Add(Array<SparseMatrix *> & Ai);
796DenseMatrix *
OuterProduct(
const DenseMatrix &A,
const DenseMatrix &B);
799SparseMatrix *
OuterProduct(
const DenseMatrix &A,
const SparseMatrix &B);
802SparseMatrix *
OuterProduct(
const SparseMatrix &A,
const DenseMatrix &B);
805SparseMatrix *
OuterProduct(
const SparseMatrix &A,
const SparseMatrix &B);
817 for (
int i = 0; i <
width; i++)
822 for (
RowNode *node_p =
Rows[row]; node_p != NULL; node_p = node_p->
Prev)
832 for (
int i = 0; i <
width; i++)
837 for (
int j =
I[row], end =
I[row+1]; j < end; j++)
850 node_p = node_p->
Prev)
871#ifdef MFEM_USE_MEMALLOC
881 return node_p->
Value;
886 MFEM_VERIFY(j != -1,
"Entry for column " << col <<
" is not allocated.");
896 return (node_p == NULL) ? 0.0 : node_p->
Value;
901 return (j == -1) ? 0.0 :
A[j];
911 for (node_p =
Rows[row]; 1; node_p = node_p->
Prev)
915#ifdef MFEM_USE_MEMALLOC
926 else if (node_p->
Column == col)
931 return node_p->
Value;
935 int *Ip =
I+row, *Jp =
J;
936 for (
int k = Ip[0], end = Ip[1]; k < end; k++)
943 MFEM_ABORT(
"Could not find entry for row = " << row <<
", col = " << col);
Abstract data type for sparse matrices.
Dynamic 2D array using row-major layout.
Data type dense matrix using column-major storage.
static MemoryClass GetHostMemoryClass()
Get the current Host MemoryClass. This is the MemoryClass used by most MFEM host Memory objects.
static MemoryClass GetDeviceMemoryClass()
Get the current Device MemoryClass. This is the MemoryClass used by most MFEM device kernels to acces...
Abstract data type for matrix inverse.
Class used by MFEM to store pointers to host and/or device memory.
void SetHostPtrOwner(bool own) const
Set/clear the ownership flag for the host pointer. Ownership indicates whether the pointer will be de...
int Capacity() const
Return the size of the allocated memory.
bool OwnsHostPtr() const
Return true if the host pointer is owned. Ownership indicates whether the pointer will be deleted by ...
bool Empty() const
Return true if the Memory object is empty, see Reset().
int width
Dimension of the input / number of columns in the matrix.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
DiagonalPolicy
Defines operator diagonal policy upon elimination of rows and/or columns.
@ DIAG_ONE
Set the diagonal value to one.
@ DIAG_ZERO
Set the diagonal value to zero.
Type
Enumeration defining IDs for some classes derived from Operator.
@ MFEM_SPARSEMAT
ID for class SparseMatrix.
int GetRow(const int row, Array< int > &cols, Vector &srow) const override
Extract all column indices and values from a given row.
void PrintCSR2(std::ostream &out) const
Prints a sparse matrix to stream out in CSR format.
const int * HostReadJ() const
void GetDiag(Vector &d) const
Returns the Diagonal of A.
int * ReadWriteI(bool on_dev=true)
real_t GetRowNorml1(int irow) const
For i = irow compute .
void PartMult(const Array< int > &rows, const Vector &x, Vector &y) const
bool isSorted
Are the columns sorted already.
bool ColumnsAreSorted() const
Returns whether or not the columns are sorted.
void Gauss_Seidel_back(const Vector &x, Vector &y) const
bool RowIsEmpty(const int row) const
void ClearColPtr() const
Reset the "current row" set by calling SetColPtr(). This method must be called between any two calls ...
MatrixInverse * Inverse() const override
This virtual method is not supported: it always returns NULL.
void Jacobi3(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
real_t * ReadWriteData(bool on_dev=true)
Memory< real_t > A
Array with size I[height], containing the actual entries of the sparse matrix, as indexed by the I ar...
void EliminateRowCol(int rc, const real_t sol, Vector &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Eliminate row rc and column rc and modify the rhs using sol.
void SetDataOwner(bool owna)
Set the data ownership flag (A array).
SparseMatrix & operator*=(real_t a)
void MultTranspose(const Vector &x, Vector &y) const override
Multiply a vector with the transposed matrix. y = At * x.
hipsparseDnVecDescr_t vecY_descr
void SetWidth(int width_=-1)
Change the width of a SparseMatrix.
MemAlloc< RowNode, 1024 > RowNodeAlloc
bool Empty() const
Check if the SparseMatrix is empty.
static hipsparseHandle_t handle
void GetSubMatrix(const Array< int > &rows, const Array< int > &cols, DenseMatrix &subm) const
SparseMatrix * At
Transpose of A. Owned. Used to perform MultTranspose() on devices.
void PartAddMult(const Array< int > &rows, const Vector &x, Vector &y, const real_t a=1.0) const
const int * GetJ() const
Return the array J, const version.
void GetRowSums(Vector &x) const
For all i compute .
int NumNonZeroElems() const override
Returns the number of the nonzero elements in the matrix.
bool Finalized() const
Returns whether or not CSR format has been finalized.
void EliminateBC(const Array< int > &ess_dofs, DiagonalPolicy diag_policy)
Eliminate essential (Dirichlet) boundary conditions.
void PrintInfo(std::ostream &out) const
Print various sparse matrix statistics.
void MoveDiagonalFirst()
Move the diagonal entry to the first position in each row, preserving the order of the rest of the co...
void Add(const int i, const int j, const real_t val)
void EliminateRowColMultipleRHS(int rc, const Vector &sol, DenseMatrix &rhs, DiagonalPolicy dpolicy=DIAG_ONE)
Similar to EliminateRowCol(int, const double, Vector &, DiagonalPolicy), but multiple values for elim...
void Jacobi(const Vector &b, const Vector &x0, Vector &x1, real_t sc, bool use_abs_diag=false) const
void SetColPtr(const int row) const
Initialize the SparseMatrix for fast access to the entries of the given row which becomes the "curren...
void SetDiagIdentity()
If a row contains only one diag entry of zero, set it to 1.
MFEM_DEPRECATED void UseCuSparse(bool useCuSparse_=true)
Deprecated equivalent of UseGPUSparse().
const Memory< int > & GetMemoryJ() const
void AbsMultTranspose(const Vector &x, Vector &y) const
y = |At| * x, using entry-wise absolute values of the transpose of matrix A
void Gauss_Seidel_forw(const Vector &x, Vector &y) const
Gauss-Seidel forward and backward iterations over a vector x.
void MakeRef(const SparseMatrix &master)
Clear the contents of the SparseMatrix and make it a reference to master.
int CheckFinite() const
Count the number of entries that are NOT finite, i.e. Inf or Nan.
void UseGPUSparse(bool useGPUSparse_=true)
Runtime option to use cuSPARSE or hipSPARSE. Only valid when using a CUDA or HIP backend.
int * WriteJ(bool on_dev=true)
Memory< int > J
Array with size I[height], containing the column indices for all matrix entries, as indexed by the I ...
cusparseMatDescr_t matA_descr
void EliminateCol(int col, DiagonalPolicy dpolicy=DIAG_ZERO)
Eliminates the column col from the matrix.
void LoseData()
Lose the ownership of the graph (I, J) and data (A) arrays.
void EnsureMultTranspose() const
Ensures that the matrix is capable of performing MultTranspose(), AddMultTranspose(),...
void Print(std::ostream &out=mfem::out, int width_=4) const override
Prints matrix to stream out.
void _Add_(const int col, const real_t a)
Add a value to an entry in the "current row". See SetColPtr().
SparseMatrix()
Create an empty SparseMatrix.
cusparseSpMatDescr_t matA_descr
real_t & Elem(int i, int j) override
Returns reference to a_{ij}.
virtual void PrintMathematica(std::ostream &out=mfem::out) const
Prints matrix as a SparseArray for importing into Mathematica.
void AbsMult(const Vector &x, Vector &y) const
y = |A| * x, using entry-wise absolute values of matrix A
void ClearGPUSparse()
Clear the cuSPARSE/hipSPARSE descriptors. This must be called after releasing the device memory of A.
void Threshold(real_t tol, bool fix_empty_rows=false)
Remove entries smaller in absolute value than a given tolerance tol. If fix_empty_rows is true,...
Memory< int > & GetMemoryI()
void SetSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void PrintMatlab(std::ostream &out=mfem::out) const override
Prints matrix in matlab format.
const int * ReadI(bool on_dev=true) const
real_t InnerProduct(const Vector &x, const Vector &y) const
Compute y^t A x.
Memory< int > I
Array with size (height+1) containing the row offsets.
real_t & operator()(int i, int j)
Returns reference to A[i][j].
const real_t * HostReadData() const
void PrintMM(std::ostream &out=mfem::out) const
Prints matrix in Matrix Market sparse format.
void ScaleColumns(const Vector &sr)
this = this * diag(sr);
void _Add_(const int row, const int col, const real_t a)
void PrintCSR(std::ostream &out) const
Prints matrix to stream out in hypre_CSRMatrix format.
void Swap(SparseMatrix &other)
MFEM_DEPRECATED void ClearCuSparse()
Deprecated equivalent of ClearGPUSparse().
const Memory< real_t > & GetMemoryData() const
int * ReadWriteJ(bool on_dev=true)
cusparseDnVecDescr_t vecY_descr
void BooleanMultTranspose(const Array< int > &x, Array< int > &y) const
y = At * x, treating all entries as booleans (zero=false, nonzero=true).
RowNode ** Rows
Array of linked lists, one for every row. This array represents the linked list (LIL) storage format.
void SetGraphOwner(bool ownij)
Set the graph ownership flag (I and J arrays).
real_t & SearchRow(const int col)
Perform a fast search for an entry in the "current row". See SetColPtr().
int * WriteI(bool on_dev=true)
const Memory< int > & GetMemoryI() const
static int SparseMatrixCount
void AddSubMatrix(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void _Set_(const int row, const int col, const real_t a)
void Symmetrize()
(*this) = 1/2 ((*this) + (*this)^t)
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
void AddMultTranspose(const Vector &x, Vector &y, const real_t a=1.0) const override
y += At * x (default) or y += a * At * x
int * GetRowColumns(const int row)
Return a pointer to the column indices in a row.
real_t GetJacobiScaling() const
Determine appropriate scaling for Jacobi iteration.
int RowSize(const int i) const
Returns the number of elements in row i.
real_t * HostReadWriteData()
int CountSmallElems(real_t tol) const
Count the number of entries with |a_ij| <= tol.
const int * ReadJ(bool on_dev=true) const
int MaxRowSize() const
Returns the maximum number of elements among all rows.
real_t _Get_(const int col) const
Read the value of an entry in the "current row". See SetColPtr().
void AddMult(const Vector &x, Vector &y, const real_t a=1.0) const override
y += A * x (default) or y += a * A * x
virtual ~SparseMatrix()
Destroys sparse matrix.
SparseMatrix & operator=(const SparseMatrix &rhs)
Assignment operator: deep copy.
void Mult(const Vector &x, Vector &y) const override
Matrix vector multiplication.
const real_t * GetData() const
Return the element data, i.e. the array A, const version.
void OverrideSize(int height_, int width_)
Sets the height and width of the matrix.
static cusparseHandle_t handle
void EliminateRow(int row, const real_t sol, Vector &rhs)
Eliminates a column from the transpose matrix.
bool OwnsGraph() const
Get the graph ownership flag (I and J arrays).
hipsparseSpMatDescr_t matA_descr
void GetBlocks(Array2D< SparseMatrix * > &blocks) const
SparseMatrix & operator+=(const SparseMatrix &B)
Add the sparse matrix 'B' to '*this'. This operation will cause an error if '*this' is finalized and ...
void Clear()
Clear the contents of the SparseMatrix.
void ResetTranspose() const
void EliminateRowColDiag(int rc, real_t value)
Perform elimination and set the diagonal entry to the given value.
real_t * GetRowEntries(const int row)
Return a pointer to the entries in a row.
void ScaleRows(const Vector &sl)
this = diag(sl) * this;
void SetRow(const int row, const Array< int > &cols, const Vector &srow)
bool OwnsData() const
Get the data ownership flag (A array).
const int * GetI() const
Return the array I, const version.
void DiagScale(const Vector &b, Vector &x, real_t sc=1.0, bool use_abs_diag=false) const
x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
real_t * GetData()
Return the element data, i.e. the array A.
void SetSubMatrixTranspose(const Array< int > &rows, const Array< int > &cols, const DenseMatrix &subm, int skip_zeros=1)
void BuildTranspose() const
Build and store internally the transpose of this matrix which will be used in the methods AddMultTran...
MemoryClass GetMemoryClass() const override
Return the MemoryClass preferred by the Operator.
void SortColumnIndices()
Sort the column indices corresponding to each row.
void Finalize(int skip_zeros=1) override
Finalize the matrix initialization, switching the storage format from LIL to CSR.
real_t IsSymmetric() const
Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix.
DenseMatrix * ToDenseMatrix() const
Produces a DenseMatrix from a SparseMatrix.
void _Set_(const int col, const real_t a)
Set an entry in the "current row". See SetColPtr().
void EliminateZeroRows(const real_t threshold=1e-12) override
If a row contains only zeros, set its diagonal to 1.
void EliminateCols(const Array< int > &cols, const Vector *x=NULL, Vector *b=NULL)
Eliminate all columns i for which cols[i] != 0.
int * GetJ()
Return the array J.
Memory< int > & GetMemoryJ()
Memory< real_t > & GetMemoryData()
cusparseDnVecDescr_t vecX_descr
const real_t * ReadData(bool on_dev=true) const
void ScaleRow(const int row, const real_t scale)
real_t * WriteData(bool on_dev=true)
hipsparseDnVecDescr_t vecX_descr
void AddRow(const int row, const Array< int > &cols, const Vector &srow)
void Jacobi2(const Vector &b, const Vector &x0, Vector &x1, real_t sc=1.0) const
int * GetI()
Return the array I.
int ActualWidth() const
Returns the actual Width of the matrix.
const int * HostReadI() const
void Set(const int i, const int j, const real_t val)
int Size() const
For backward compatibility, define Size() to be synonym of Height().
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
DenseMatrix * OuterProduct(const DenseMatrix &A, const DenseMatrix &B)
Produces a block matrix with blocks A_{ij}*B.
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
void Swap< SparseMatrix >(SparseMatrix &a, SparseMatrix &b)
Specialization of the template function Swap<> for class SparseMatrix.
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
T * Write(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for write access to mem with the mfem::Device's DeviceMemoryClass, if on_dev = true,...
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
MemoryClass
Memory classes identify sets of memory types.
T * ReadWrite(Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read+write access to mem with the mfem::Device's DeviceMemoryClass,...
void RAP(const DenseMatrix &A, const DenseMatrix &P, DenseMatrix &RAP)
SparseMatrix * MultAbstractSparseMatrix(const AbstractSparseMatrix &A, const AbstractSparseMatrix &B)
Matrix product of sparse matrices. A and B do not need to be CSR matrices.
SparseMatrix * TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A, int useActualWidth)
Transpose of a sparse matrix. A does not need to be a CSR matrix.
SparseMatrix * Mult_AtDA(const SparseMatrix &A, const Vector &D, SparseMatrix *OAtDA)
Matrix multiplication A^t D A. All matrices must be finalized.
MemoryType
Memory types supported by MFEM.
std::function< real_t(const Vector &)> f(real_t mass_coeff)
SparseMatrix * TransposeMult(const SparseMatrix &A, const SparseMatrix &B)
C = A^T B.
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
void SparseMatrixFunction(SparseMatrix &S, real_t(*f)(real_t))
Applies f() to each element of the matrix (after it is finalized).