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);
  int          (*zero)(SUNMatrix);
  int          (*copy)(SUNMatrix, SUNMatrix);
  int          (*scaleadd)(realtype, SUNMatrix, SUNMatrix);
  int          (*scaleaddi)(realtype, SUNMatrix);
  int          (*matvecsetup)(SUNMatrix);
  int          (*matvec)(SUNMatrix, N_Vector, N_Vector);
  int          (*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:

int SUNMatZero(SUNMatrix A)
  return((int) 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.

int SUNMatCopyOps(SUNMatrix A, SUNMatrix B)

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

  • A – the matrix to copy operations from.

  • B – the matrix to copy operations to.

Return value:

If successful, this function returns 0. If either of the inputs are NULL or the ops structure of either input is NULL, then is function returns a non-zero value.

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.

  • 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


Band \(M \times M\) matrix


CUDA sparse CSR matrix


User-provided custom matrix


Dense \(M \times N\) matrix


SUNMatrix wraper for Ginkgo matrices


Dense \(M \times N\) matrix


oneMKL dense \(M \times N\) matrix


SUNMatrix wrapper for SuperLU_DIST SuperMatrix


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


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.


B = SUNMatClone(A);
void SUNMatDestroy(SUNMatrix A)

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


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

Returns the storage requirements for the matrix A. lrw contains the number of realtype 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.


retval = SUNMatSpace(A, &lrw, &liw);
int SUNMatZero(SUNMatrix A)

Zeros all entries of the SUNMatrix A. The return value is an integer flag denoting success/failure of the operation:

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


retval = SUNMatZero(A);
int SUNMatCopy(SUNMatrix A, SUNMatrix B)

Performs the operation B gets A for all entries of the matrices A and B. The return value is an integer flag denoting success/failure of the operation:

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


retval = SUNMatCopy(A,B);
int SUNMatScaleAdd(realtype c, SUNMatrix A, SUNMatrix B)

Performs the operation A gets cA + B. The return value is an integer flag denoting success/failure of the operation:

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


retval = SUNMatScaleAdd(c, A, B);
int SUNMatScaleAddI(realtype c, SUNMatrix A)

Performs the operation A gets cA + I. The return value is an integer flag denoting success/failure of the operation:

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


retval = SUNMatScaleAddI(c, A);
int SUNMatMatvecSetup(SUNMatrix A)

Performs any setup necessary to perform a matrix-vector product. The return value is an integer flag denoting 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.


retval = SUNMatMatvecSetup(A);
int 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 is an integer flag denoting success/failure of the operation:

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


retval = SUNMatMatvec(A, x, y);

10.2.1. SUNMatrix return codes

The functions provided to SUNMatrix modules within the SUNDIALS-provided SUNMatrix implementations utilize a common set of return codes, listed below. These adhere to a common pattern: 0 indicates success, a negative value indicates a failure. Aside from this pattern, the actual values of each error code are primarily to provide additional information to the user in case of a SUNMatrix failure.

  • SUNMAT_SUCCESS (0) – successful call

  • SUNMAT_ILL_INPUT (-1) – an illegal input has been provided to the function

  • SUNMAT_MEM_FAIL (-2) – failed memory access or allocation

  • SUNMAT_OPERATION_FAIL (-3) – a SUNMatrix operation returned nonzero

  • SUNMAT_MATVEC_SETUP_REQUIRED (-4) – the SUNMatMatvecSetup() routine needs to be called prior to calling SUNMatMatvec()