Working with sparse data in RandBLAS

Sparse matrix data structures

The common interface for our sparse matrix types

enum class RandBLAS::sparse_data::IndexBase : int

Indicates whether the rows and/or columns of a sparse matrix are enumerated (that is, indexed) starting from zero or starting from one. The majority of RandBLAS’ sparse matrix functionality requires zero-based indexing.

Values:

enumerator Zero

enumerator One

template<typename SpMat> concept SparseMatrix

#include <base.hh>

An object $\texttt{M}$ of type $\texttt{SpMat}$ has the following attributes.

	type	description
$\texttt{M.n_rows}$	$\texttt{const int64_t}$	number of rows
$\texttt{M.n_cols}$	$\texttt{const int64_t}$	number of columns
$\texttt{M.nnz}$	$\texttt{int64_t}$	number of structural nonzeros
$\texttt{M.vals}$	$\texttt{SpMat::scalar_t *}$	pointer to values of structural nonzeros
$\texttt{M.own_memory}$	$\texttt{bool}$	A flag indicating if memory attached to $\texttt{M}$ should be deallocated when $\texttt{M}$ is deleted. This flag is set automatically based on the type of constructor used for $\texttt{M}.$

Memory-owning constructors

$\texttt{SpMat}$ must have a constructor for an empty matrix of given dimensions. Conformant implementations of this constructor look like the following.
SpMat(int64_t n_rows, int64_t n_cols) 
 : n_rows(n_rows), n_cols(n_cols), nnz(0), vals(nullptr), own_memory(true) {
     // class-specific code ...
 };
If we construct $\texttt{SpMat M(m, n)},$ then we can’t store data in $\texttt{M}$ until a function call of the form $\texttt{M.reserve(nnz)}.$ Here’s an outline of a conformant implementation of this function.
void reserve(int64_t nnz) {
    assert this->own_memory;
    this->nnz = nnz;
    this->vals = new SpMat::scalar_t[nnz];
    // ... class-specific code ...
}
The destructor of $\texttt{M}$ is responsible for deallocating $\texttt{M.vals}$ and other attached data. A conformant implemnentation of the destructor will look like the following.
~SpMat() {
    if (own_memory) {
        if (vals != nullptr) delete [] vals;
        // ... class-specific code ...
    }
}        

Instance methods

In addition to $\texttt{SpMat::reserve}$, this concept requires two more instance methods. Let $\texttt{A}$ denote an $\texttt{SpMat}.$

Calling $\texttt{A.deepcopy()}$ returns an $\texttt{SpMat}$ that’s mathematically equivalent to $\texttt{A}$ and that owns its attached memory.

Calling $\texttt{A.transpose()}$ returns an instance of a (const) type that conforms to the SparseMatrix concept and that is mathematically equivalent to the transpose of $\texttt{A}.$

View constructors

This concept doesn’t place requirements on constructors for sparse matrix views of existing data. However, all of RandBLAS’ sparse matrix classes offer such constructors. See individual classes’ documentation for details.

Passing to and returning from functions

As a consequence of our requirements on $\texttt{SpMat::deepcopy}$, $\texttt{SpMat}$ must have a C++ move constructor. The move constructors for RandBLAS’ SparseMatrix types do not appear in our web documentation since move constructors should not be called by user code.

$\texttt{SpMat}$ does not necessarily have a C++ copy constructor. In fact, none of RandBLAS’ SparseMatrix classes have a C++ copy constructor. If $\texttt{SpMat}$ has no copy constructor then instances of that type can only be passed by reference.

COOMatrix

enum class RandBLAS::sparse_data::NonzeroSort : char

Indicates whether the (vals, rows, cols) data of a COO-format sparse matrix are known to be sorted in CSC order, CSR order, or neither of those orders.

Values:

enumerator CSC

enumerator CSR

enumerator None

template<typename T, SignedInteger sint_t = int64_t> struct COOMatrix

Let ${\mathbf{A}}$ denote a sparse matrix with ${\texttt{nnz}}$ structural nonzeros. Its COO representation consists of declared dimensions, ${\texttt{n_rows}}$ and ${\texttt{n_cols}}$ , as well as a triplet of arrays ${(\texttt{vals},\texttt{rows},\texttt{cols})}$ where

\[\mathbf{A}_{\texttt{rows}[\ell],\texttt{cols}[\ell]} = \texttt{vals}[\ell] \quad\text{for all}\quad \ell \in \{0,\ldots,\texttt{nnz}-1\}.\]

This type conforms to the SparseMatrix concept.

Public Types

using scalar_t = T : Real scalar type used for structural nonzeros in this matrix.

using index_t = sint_t : Signed integer type used in the rows and cols array members.

Public Functions

inline COOMatrix(int64_t n_rows, int64_t n_cols)

Standard constructor. Initializes n_rows and n_cols at the provided values. The vals, rows, and cols members are set to null pointers; nnz is set to Zero, index_base is set to zero, and COOMatrix::own_memory is set to true.

This constructor is intended for use with COOMatrix::reserve(int64_t nnz).

inline COOMatrix(int64_t n_rows, int64_t n_cols, int64_t nnz, T *vals, sint_t *rows, sint_t *cols, bool compute_sort_type = true, IndexBase index_base = IndexBase::Zero)

Expert constructor. Arguments passed to this function are used to initialize members of the same names; COOMatrix::own_memory is set to false.

If compute_sort_type is true, then COOMatrix::sort will be computed by inspecting the contents of (rows, cols). Otherwise, COOMatrix::sort is set to None.

inline void reserve(int64_t arg_nnz)

This function requires that own_memory is true, that arg_nnz > 0, and that vals, rows, and cols are null. If any of these conditions are not met then this function will raise an error.

If no error is raised then this function redirects vals, rows, and cols to new arrays of length arg_nnz, and nnz is set to arg_nnz.

inline void sort_arrays(NonzeroSort s)

Sort the (vals, rows, cols) underlying this matrix to facilitate fast conversion to CSR or CSC formats.

This function has no effect if s == COOMatrix::sort or s == NonzeroSort::None.

template<SignedInteger index_t> inline void symperm_inplace(const index_t *perm, bool preserve_sort = true)

This function requires that n_rows = n_cols and and that ${\texttt{index_base}}$ is Zero.

Define ${n}$ := n_rows. The buffer ${\texttt{perm}}$ is a permutation of ${\{0, 1,\ldots, n - 1\}}$ that defines an ${n \times n}$ permutation matrix ${\mathbf{P} = \mathbf{I}_n(\texttt{perm},:).}$ This function overwrites the current sparse matrix ${\mathbf{A}}$ by ${\mathbf{A} := \mathbf{P} \mathbf{A} \mathbf{P}^T.}$ If ${\texttt{preserve_sort}}$ is true then we ensure that ${\mathbf{A}.\texttt{sort}}$ has the same value on exit as it had on entry. Otherwise, we overwrite ${\mathbf{A}.\texttt{sort}}$ with None.

inline COOMatrix<T, sint_t> deepcopy() const: Return a memory-owning copy of this COOMatrix.

inline CSRMatrix<T, sint_t> as_owning_csr() const

Return a memory-owning CSRMatrix repesentation of this COOMatrix.

If sort != NonzeroSort::CSR, then we create a temporary deep copy of this matrix on the way to constructing the CSR representation.

inline CSCMatrix<T, sint_t> as_owning_csc() const

Return a memory-owning CSCMatrix representation of this COOMatrix.

If sort != NonzeroSort::CSC, then we create a temporary deep copy of this matrix on the way to constructing the CSC representation.

inline const COOMatrix<T, sint_t> transpose() const: Return a const view of the transpose of this COOMatrix.

Public Members

const int64_t n_rows: The number of rows in this sparse matrix.

const int64_t n_cols: The number of columns in this sparse matrix.

bool own_memory

If true, then RandBLAS has permission to allocate and attach memory to the reference members of this matrix (vals, rows, and cols). If true at destruction time, then delete [] will be called on each non-null reference member of this matrix.

RandBLAS only writes to this member at construction time.

int64_t nnz: The number of structral nonzeros in this matrix.

IndexBase index_base: A flag to indicate whether (rows, cols) are interpreted with zero-based or one-based indexing.

T *vals

Reference to an array that holds values for the structural nonzeros of this matrix.

If non-null, this must have length at least nnz.

sint_t *rows

Reference to an array that holds row indices for the structural nonzeros of this matrix.

If non-null, this must have length at least nnz.

sint_t *cols

Reference to an array that holds column indicies for the structural nonzeros of this matrix.

If non-null, this must have length at least nnz.

NonzeroSort sort: A flag to indicate if the data in (vals, rows, cols) is sorted in a CSC-like order, a CSR-like order, or neither order.

CSRMatrix

template<typename T, SignedInteger sint_t = int64_t> struct CSRMatrix

Let ${\mathbf{A}}$ denote a sparse matrix with ${\texttt{nnz}}$ structural nonzeros. Its CSR representation consists of declared dimensions, ${\texttt{n_rows}}$ and ${\texttt{n_cols}}$ , and a triplet of arrays ${(\texttt{vals},\texttt{rowptr},\texttt{colidxs}).}$

The ${\texttt{i}^{\text{th}}}$ row of ${\mathbf{A}}$ has ${\texttt{rowptr[i+1] - rowptr[i]}}$ structural nonzeros. The ${\texttt{k}^{\text{th}}}$ structural nonzero in row ${\texttt{i}}$ appears in column ${\texttt{colidxs[rowptr[i] + k]}}$ and is equal to ${\texttt{vals[rowptr[i] + k]}.}$

This type conforms to the SparseMatrix concept.

Public Types

using scalar_t = T : Real scalar type used for structural nonzeros in this matrix.

using index_t = sint_t : Signed integer type used in the rowptr and colidxs array members.

Public Functions

inline CSRMatrix(int64_t n_rows, int64_t n_cols)

Standard constructor. Initializes n_rows and n_cols at the provided values. The vals, rowptr, and colidxs members are set to null pointers; nnz is set to zero, index_base is set to Zero, and CSRMatrix::own_memory is set to true.

This constructor is intended for use with CSRMatrix::reserve(int64_t arg_nnz).

inline CSRMatrix(int64_t n_rows, int64_t n_cols, int64_t nnz, T *vals, sint_t *rowptr, sint_t *colidxs, IndexBase index_base = IndexBase::Zero): Expert constructor. Arguments passed to this function are used to initialize members of the same names; CSRMatrix::own_memory is set to false.

inline void reserve(int64_t arg_nnz)

This function requires that arg_nnz > 0, that own_memory is true, that colidxs and vals are null. If any of these conditions are not met then this function will raise an error.

Special logic applies rowptr because its documented length requirement is equal to the const expression (n_rows + 1).

If rowptr is non-null then it is left unchanged, and it is presumed to point to an array of length at least n_rows + 1.
If rowptr is null, then it will be redirected to a new array of type sint_t and length n_rows + 1.

From there, the reference members colidxs and vals are redirected to new arrays length arg_nnz, and nnz is updated to arg_nnz.

inline CSRMatrix<T, sint_t> deepcopy() const: Return a memory-owning copy of this CSRMatrix.

inline COOMatrix<T, sint_t> as_owning_coo() const: Return a memory-owning COOMatrix representation of this CSRMatrix.

inline const CSCMatrix<T, sint_t> transpose() const: Return a const CSCMatrix view of the transpose of this CSRMatrix.

Public Members

const int64_t n_rows: The number of rows in this sparse matrix.

const int64_t n_cols: The number of columns in this sparse matrix.

bool own_memory

If true, then RandBLAS has permission to allocate and attach memory to the reference members of this matrix (vals, rowptr, colidxs). If true at destruction time, then delete [] will be called on each non-null reference member of this matrix.

RandBLAS only writes to this member at construction time.

int64_t nnz: The number of structral nonzeros in this matrix.

IndexBase index_base: A flag to indicate whether colidxs is interpreted with zero-based or one-based indexing.

T *vals

Reference to an array that holds values for the structural nonzeros of this matrix.

If non-null, this must have length at least nnz.

sint_t *rowptr

Reference to a pointer offset array for the CSR format.

If non-null, then must have length at least ${\texttt{n_rows + 1}}$ .

sint_t *colidxs

Reference to a column index array in the CSR format, interpreted in ${\texttt{index_base}}$ .

If non-null, then must have length at least nnz.

CSCMatrix

template<typename T, SignedInteger sint_t = int64_t> struct CSCMatrix

Let ${\mathbf{A}}$ denote a sparse matrix with ${\texttt{nnz}}$ structural nonzeros. Its CSC representation consists of declared dimensions, ${\texttt{n_rows}}$ and ${\texttt{n_cols}}$ , and a triplet of arrays ${(\texttt{vals},\texttt{rowidxs},\texttt{colptr}).}$

The ${\texttt{j}^{\text{th}}}$ column of ${\mathbf{A}}$ has ${\texttt{colptr[j+1] - colptr[j]}}$ structural nonzeros. The ${\texttt{k}^{\text{th}}}$ structural nonzero in column ${\texttt{j}}$ appears in row ${\texttt{rowidxs[colptr[j] + k]}}$ and is equal to ${\texttt{vals[colptr[j] + k]}.}$

This type conforms to the SparseMatrix concept.

Public Types

using scalar_t = T : Real scalar type used for structural nonzeros in this matrix.

using index_t = sint_t : Signed integer type used in the rowidxs and colptr array members.

Public Functions

inline CSCMatrix(int64_t n_rows, int64_t n_cols)

Standard constructor. Initializes n_rows and n_cols at the provided values. The vals, rowidxs, and colptr members are null-initialized; ${\texttt{nnz}}$ is set to zero, ${\texttt{index_base}}$ is set to Zero, and CSCMatrix::own_memory is set to true.

This constructor is intended for use with COOMatrix::reserve(int64_t nnz).

inline CSCMatrix(int64_t n_rows, int64_t n_cols, int64_t nnz, T *vals, sint_t *rowidxs, sint_t *colptr, IndexBase index_base = IndexBase::Zero): Expert constructor. Arguments passed to this function are used to initialize members of the same names; CSCMatrix::own_memory is set to false.

inline void reserve(int64_t arg_nnz)

This function requires that arg_nnz > 0, that own_memory is true, that rowidxs and vals are null. If any of these conditions are not met then this function will raise an error.

Special logic applies colptr because its documented length requirement is equal to the const expression (n_cols + 1).

If colptr is non-null then it is left unchanged, and it is presumed to point to an array of length at least n_cols + 1.
If colptr is null, then it will be redirected to a new array of type sint_t and length n_cols + 1.

From there, the reference members rowidxs and vals are redirected to new arrays length arg_nnz, and nnz is updated to arg_nnz.

inline CSCMatrix<T, sint_t> deepcopy() const: Return a memory-owning copy of this CSCMatrix.

inline COOMatrix<T, sint_t> as_owning_coo() const: Return a memory-owning COOMatrix representation of this CSCMatrix.

inline const CSRMatrix<T, sint_t> transpose() const: Return a const CSRMatrix view of the transpose of this CSCMatrix.

Public Members

const int64_t n_rows: The number of rows in this sparse matrix.

const int64_t n_cols: The number of columns in this sparse matrix.

bool own_memory

If true, then RandBLAS has permission to allocate and attach memory to the reference members of this matrix (vals, rowidxs, colptr). If true at destruction time, then delete [] will be called on each non-null reference member of this matrix.

RandBLAS only writes to this member at construction time.

int64_t nnz: The number of structral nonzeros in this matrix.

IndexBase index_base: A flag to indicate whether rowidxs is interpreted with zero-based or one-based indexing.

T *vals

Reference to an array that holds values for the structural nonzeros of this matrix.

If non-null, this must have length at least nnz.

sint_t *rowidxs

Reference to a row index array in the CSC format, interpreted in ${\texttt{index_base}}$ .

If non-null, then must have length at least nnz.

sint_t *colptr

Reference to a pointer offset array for the CSC format.

If non-null, then must have length at least ${\texttt{n_cols + 1}}$ .

Operations with sparse matrices

Sketching

$\mathbf{B} = \alpha \cdot \operatorname{op}(\operatorname{submat}(\mathbf{S}))\cdot \operatorname{op}(\mathbf{A}) + \beta \cdot \mathbf{B}$

template<SparseMatrix SpMat, typename DenseSkOp, typename T = DenseSkOp::scalar_t> inline void RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opS, blas::Op opA, int64_t d, int64_t n, int64_t m, T alpha, const DenseSkOp &S, int64_t ro_s, int64_t co_s, const SpMat &A, T beta, T *B, int64_t ldb)

Sketch from the left in an SpMM-like operation

\[\operatorname{mat}(B) = \alpha \cdot \underbrace{\operatorname{op}(\operatorname{submat}(\mathbf{S}))}_{d \times m} \cdot \underbrace{\operatorname{op}(\mathbf{A})}_{m \times n} + \beta \cdot \underbrace{\operatorname{mat}(B)}_{d \times n}, \tag{$\star$}\]

where $\alpha$ and $\beta$ are real scalars, $\operatorname{op}(\mathbf{X})$ either returns a matrix $\mathbf{X}$ or its transpose, $\mathbf{A}$ is a sparse matrix, and $\mathbf{S}$ is a dense sketching operator.

$\mathbf{B} = \alpha \cdot \operatorname{op}(\mathbf{A})\cdot \operatorname{op}(\operatorname{submat}(\mathbf{S})) + \beta \cdot \mathbf{B}$

template<SparseMatrix SpMat, typename DenseSkOp, typename T = DenseSkOp::scalar_t> inline void RandBLAS::sketch_sparse(blas::Layout layout, blas::Op opA, blas::Op opS, int64_t m, int64_t d, int64_t n, T alpha, const SpMat &A, const DenseSkOp &S, int64_t ro_s, int64_t co_s, T beta, T *B, int64_t ldb)

Sketch from the right in an SpMM-like operation

\[\operatorname{mat}(B) = \alpha \cdot \underbrace{\operatorname{op}(\mathbf{A})}_{m \times n} \cdot \underbrace{\operatorname{op}(\operatorname{submat}(\mathbf{S}))}_{n \times d} + \beta \cdot \underbrace{\operatorname{mat}(B)}_{m \times d}, \tag{$\star$}\]

where $\alpha$ and $\beta$ are real scalars, $\operatorname{op}(\mathbf{X})$ either returns a matrix $\mathbf{X}$ or its transpose, $\mathbf{A}$ is a sparse matrix, and $\mathbf{S}$ is a dense sketching operator.

Deterministic operations

$\mathbf{C} = \alpha \cdot \operatorname{op}(\mathbf{A})\cdot \operatorname{op}(\mathbf{B}) + \beta \cdot \mathbf{C},$ with sparse $\mathbf{A}$

template<SparseMatrix SpMat, typename T = SpMat::scalar_t> inline void RandBLAS::spmm(blas::Layout layout, blas::Op opA, blas::Op opB, int64_t m, int64_t n, int64_t k, T alpha, const SpMat &A, const T *B, int64_t ldb, T beta, T *C, int64_t ldc)

Multiply a dense matrix on the left with a sparse matrix:

\[\operatorname{mat}(C) = \alpha \cdot \underbrace{\operatorname{op}(\mathbf{A})}_{m \times k} \cdot \underbrace{\operatorname{op}(\operatorname{mat}(B))}_{k \times n} + \beta \cdot \underbrace{\operatorname{mat}(C)}_{m \times n}, \tag{$\star$}\]

where $\alpha$ and $\beta$ are real scalars, $\operatorname{op}(\mathbf{X})$ either returns a matrix $\mathbf{X}$ or its transpose, and $\mathbf{A}$ is sparse.

$\mathbf{C} = \alpha \cdot \operatorname{op}(\mathbf{A})\cdot \operatorname{op}(\mathbf{B}) + \beta \cdot \mathbf{C},$ with sparse $\mathbf{B}$

template<SparseMatrix SpMat, typename T = SpMat::scalar_t> inline void RandBLAS::spmm(blas::Layout layout, blas::Op opA, blas::Op opB, int64_t m, int64_t n, int64_t k, T alpha, const T *A, int64_t lda, const SpMat &B, T beta, T *C, int64_t ldc)

Multiply a dense matrix on the right with a sparse matrix:

\[\operatorname{mat}(C) = \alpha \cdot \underbrace{\operatorname{op}(\operatorname{mat}(A))}_{m \times k} \cdot \underbrace{\operatorname{op}(\mathbf{B})}_{k \times n} + \beta \cdot \underbrace{\operatorname{mat}(C)}_{m \times n}, \tag{$\star$}\]

where $\alpha$ and $\beta$ are real scalars, $\operatorname{op}(\mathbf{X})$ either returns a matrix $\mathbf{X}$ or its transpose, and $\mathbf{B}$ is sparse.

$\mathbf{B} = \alpha \cdot \operatorname{op}(\mathbf{A})^{-1} \cdot \mathbf{B},$ with sparse triangular $\mathbf{A}$

template<SparseMatrix SpMat, typename T = SpMat::scalar_t> void RandBLAS::sparse_data::trsm(blas::Layout layout, blas::Op opA, T alpha, const SpMat &A, blas::Uplo uplo, blas::Diag diag, int64_t n, T *B, int64_t ldb, int validation_mode = 1)

Overwrite $n$ columns of a matrix $\operatorname{mat}(B)$ with the results of a scaled triangular solve

\[\operatorname{mat}(B) = \alpha \cdot \underbrace{\operatorname{op}(\mathbf{A})^{-1}}_{m \times m} \cdot \underbrace{\operatorname{mat}(B)}_{m \times n}, \tag{$\star$}\]

where $\mathbf{A}$ is either the sparse matrix $A$ or a view of $A$ that replaces its diagonal with the vector of all ones, and $\operatorname{op}(\mathbf{A})$ returns either $\mathbf{A}$ or $\mathbf{A}^T.$

	type	description
\(\texttt{M.n_rows}\)	\(\texttt{const int64_t}\)	number of rows
\(\texttt{M.n_cols}\)	\(\texttt{const int64_t}\)	number of columns
\(\texttt{M.nnz}\)	\(\texttt{int64_t}\)	number of structural nonzeros
\(\texttt{M.vals}\)	\(\texttt{SpMat::scalar_t *}\)	pointer to values of structural nonzeros
\(\texttt{M.own_memory}\)	\(\texttt{bool}\)	A flag indicating if memory attached to \(\texttt{M}\) should be deallocated when \(\texttt{M}\) is deleted. This flag is set automatically based on the type of constructor used for \(\texttt{M}.\)