10.1. Description of the SUNMATRIX Modules

For problems that involve direct methods for solving linear systems, the SUNDIALS packages not only operate on generic vectors, but also on generic matrices (of type SUNMatrix), through a set of operations defined by the particular SUNMATRIX implementation. Users can provide their own specific implementation of the SUNMATRIX module, particularly in cases where they provide their own N_Vector and/or linear solver modules, and require matrices that are compatible with those implementations. The generic SUNMatrix operations are described below, and descriptions of the SUNMATRIX implementations provided with SUNDIALS follow.

The generic SUNMatrix type has been modeled after the object-oriented style of the generic N_Vector type. Specifically, a generic SUNMatrix is a pointer to a structure that has an implementation-dependent content field containing the description and actual data of the matrix, and an ops field pointing to a structure with generic matrix operations. The type SUNMatrix is defined as:

typedef struct _generic_SUNMatrix *SUNMatrix

and the generic structure is defined as

struct _generic_SUNMatrix {
    void *content;
    struct _generic_SUNMatrix_Ops *ops;
};

Here, the _generic_SUNMatrix_Ops structure is essentially a list of function pointers to the various actual matrix operations, and is defined as

struct _generic_SUNMatrix_Ops {
  SUNMatrix_ID (*getid)(SUNMatrix);
  SUNMatrix    (*clone)(SUNMatrix);
  void         (*destroy)(SUNMatrix);
  SUNErrCode   (*zero)(SUNMatrix);
  SUNErrCode   (*copy)(SUNMatrix, SUNMatrix);
  SUNErrCode   (*scaleadd)(sunrealtype, SUNMatrix, SUNMatrix);
  SUNErrCode   (*scaleaddi)(sunrealtype, SUNMatrix);
  SUNErrCode   (*matvecsetup)(SUNMatrix);
  SUNErrCode   (*matvec)(SUNMatrix, N_Vector, N_Vector);
  SUNErrCode   (*space)(SUNMatrix, long int*, long int*);
};

The generic SUNMATRIX module defines and implements the matrix operations acting on a SUNMatrix. These routines are nothing but wrappers for the matrix operations defined by a particular SUNMATRIX implementation, which are accessed through the ops field of the SUNMatrix structure. To illustrate this point we show below the implementation of a typical matrix operation from the generic SUNMATRIX module, namely SUNMatZero, which sets all values of a matrix A to zero, returning a flag denoting a successful/failed operation:

SUNErrCode SUNMatZero(SUNMatrix A)
{
  return(A->ops->zero(A));
}

§10.2 contains a complete list of all matrix operations defined by the generic SUNMATRIX module. A particular implementation of the SUNMATRIX module must:

Specify the content field of the SUNMatrix object.
Define and implement a minimal subset of the matrix operations. See the documentation for each SUNDIALS package and/or linear solver to determine which SUNMATRIX operations they require.

Note that the names of these routines should be unique to that implementation in order to permit using more than one SUNMATRIX module (each with different SUNMatrix internal data representations) in the same code.
Define and implement user-callable constructor and destructor routines to create and free a SUNMatrix with the new content field and with ops pointing to the new matrix operations.
Optionally, define and implement additional user-callable routines acting on the newly defined SUNMatrix (e.g., a routine to print the content for debugging purposes).
Optionally, provide accessor macros as needed for that particular implementation to be used to access different parts in the content field of the newly defined SUNMatrix.

To aid in the creation of custom SUNMATRIX modules the generic SUNMATRIX module provides three utility functions SUNMatNewEmpty(), SUNMatCopyOps(), and SUNMatFreeEmpty(). When used in custom SUNMATRIX constructors and clone routines these functions will ease the introduction of any new optional matrix operations to the SUNMATRIX API by ensuring only required operations need to be set and all operations are copied when cloning a matrix.

SUNMatrix SUNMatNewEmpty()

This function allocates a new generic SUNMatrix object and initializes its content pointer and the function pointers in the operations structure to NULL.

Return value:: If successful, this function returns a SUNMatrix object. If an error occurs when allocating the object, then this routine will return NULL.

SUNErrCode SUNMatCopyOps(SUNMatrix A, SUNMatrix B)

This function copies the function pointers in the ops structure of A into the ops structure of B.

Arguments:

A – the matrix to copy operations from.

B – the matrix to copy operations to.

Return value:

A SUNErrCode

void SUNMatFreeEmpty(SUNMatrix A)

This routine frees the generic SUNMatrix object, under the assumption that any implementation-specific data that was allocated within the underlying content structure has already been freed. It will additionally test whether the ops pointer is NULL, and, if it is not, it will free it as well.

Arguments:

A – the SUNMatrix object to free

Each SUNMATRIX implementation included in SUNDIALS has a unique identifier specified in enumeration and shown in Table 10.1. It is recommended that a user-supplied SUNMATRIX implementation use the SUNMATRIX_CUSTOM identifier.

Table 10.1 Identifiers associated with matrix kernels supplied with SUNDIALS
Matrix ID	Matrix type
SUNMATRIX_BAND	Band \(M \times M\) matrix
SUNMATRIX_CUSPARSE	CUDA sparse CSR matrix
SUNMATRIX_CUSTOM	User-provided custom matrix
SUNMATRIX_DENSE	Dense \(M \times N\) matrix
SUNMATRIX_GINKGO	SUNMatrix wraper for Ginkgo matrices
SUNMATRIX_MAGMADENSE	Dense \(M \times N\) matrix
SUNMATRIX_ONEMKLDENSE	oneMKL dense \(M \times N\) matrix
SUNMATRIX_SLUNRLOC	SUNMatrix wrapper for SuperLU_DIST SuperMatrix
SUNMATRIX_SPARSE	Sparse (CSR or CSC) \(M\times N\) matrix

10.2. Description of the SUNMATRIX operations

For each of the SUNMatrix operations, we give the name, usage of the function, and a description of its mathematical operations below.

SUNMatrix_ID SUNMatGetID(SUNMatrix A)

Returns the type identifier for the matrix A. It is used to determine the matrix implementation type (e.g. dense, banded, sparse,…) from the abstract SUNMatrix interface. This is used to assess compatibility with SUNDIALS-provided linear solver implementations. Returned values are given in Table 10.1

Usage:

id = SUNMatGetID(A);

SUNMatrix SUNMatClone(SUNMatrix A)

Creates a new SUNMatrix of the same type as an existing matrix A and sets the ops field. It does not copy the matrix values, but rather allocates storage for the new matrix.

Usage:

B = SUNMatClone(A);

void SUNMatDestroy(SUNMatrix A)

Destroys the SUNMatrix A and frees memory allocated for its internal data.

Usage:

SUNMatDestroy(A);

SUNErrCode SUNMatSpace(SUNMatrix A, long int *lrw, long int *liw)

Returns the storage requirements for the matrix A. lrw contains the number of sunrealtype words and liw contains the number of integer words. The return value denotes success/failure of the operation.

This function is advisory only, for use in determining a user’s total space requirements; it could be a dummy function in a user-supplied SUNMatrix module if that information is not of interest.

Usage:

retval = SUNMatSpace(A, &lrw, &liw);

SUNErrCode SUNMatZero(SUNMatrix A)

Zeros all entries of the SUNMatrix A. The return value denotes the success/failure of the operation:

\[A_{i,j} = 0, \quad i=1,\ldots,m, \; j=1,\ldots,n.\]

Usage:

retval = SUNMatZero(A);

SUNErrCode SUNMatCopy(SUNMatrix A, SUNMatrix B)

Performs the operation B gets A for all entries of the matrices A and B. The return value denotes the success/failure of the operation:

\[B_{i,j} = A_{i,j}, \quad i=1,\ldots,m, \; j=1,\ldots,n.\]

Usage:

retval = SUNMatCopy(A,B);

SUNErrCode SUNMatScaleAdd(sunrealtype c, SUNMatrix A, SUNMatrix B)

Performs the operation A gets cA + B. The return value denotes the success/failure of the operation:

\[A_{i,j} = cA_{i,j} + B_{i,j}, \quad i=1,\ldots,m, \; j=1,\ldots,n.\]

Usage:

retval = SUNMatScaleAdd(c, A, B);

SUNErrCode SUNMatScaleAddI(sunrealtype c, SUNMatrix A)

Performs the operation A gets cA + I. The return value denotes the success/failure of the operation:

\[A_{i,j} = cA_{i,j} + \delta_{i,j}, \quad i,j=1,\ldots,n.\]

Usage:

retval = SUNMatScaleAddI(c, A);

SUNErrCode SUNMatMatvecSetup(SUNMatrix A)

Performs any setup necessary to perform a matrix-vector product. The return value denotes the success/failure of the operation. It is useful for SUNMatrix implementations which need to prepare the matrix itself, or communication structures before performing the matrix-vector product.

Usage:

retval = SUNMatMatvecSetup(A);

SUNErrCode SUNMatMatvec(SUNMatrix A, N_Vector x, N_Vector y)

Performs the matrix-vector product y gets Ax. It should only be called with vectors x and y that are compatible with the matrix A – both in storage type and dimensions. The return value denotes the success/failure of the operation:

\[y_i = \sum_{j=1}^n A_{i,j} x_j, \quad i=1,\ldots,m.\]

Usage:

retval = SUNMatMatvec(A, x, y);