11.1. The SUNLinearSolver API

The SUNLinSol API defines several linear solver operations that enable SUNDIALS packages to utilize this API. These functions can be divided into three categories. The first are the core linear solver functions. The second consist of “set” routines to supply the linear solver with functions provided by the SUNDIALS packages and to modify solver parameters. The final group consists of “get” routines for retrieving linear solver statistics. All of these functions are defined in the header file sundials/sundials_linearsolver.h.

11.1.1. SUNLinearSolver core functions

The core linear solver functions consist of two required functions: SUNLinSolGetType() returns the linear solver type, and SUNLinSolSolve() solves the linear system \(Ax=b\).

The remaining optional functions return the solver ID (SUNLinSolGetID()), initialize the linear solver object once all solver-specific options have been set (SUNLinSolInitialize()), set up the linear solver object to utilize an updated matrix \(A\) (SUNLinSolSetup()), and destroy a linear solver object (SUNLinSolFree()).

SUNLinearSolver_Type SUNLinSolGetType(SUNLinearSolver LS)

Returns the type identifier for the linear solver LS.

Return value:

• SUNLINEARSOLVER_DIRECT (0) – the SUNLinSol module requires a matrix, and computes an “exact” solution to the linear system defined by that matrix.

• SUNLINEARSOLVER_ITERATIVE (1) – the SUNLinSol module does not require a matrix (though one may be provided), and computes an inexact solution to the linear system using a matrix-free iterative algorithm. That is it solves the linear system defined by the package-supplied ATimes routine (see SUNLinSolSetATimes() below), even if that linear system differs from the one encoded in the matrix object (if one is provided). As the solver computes the solution only inexactly (or may diverge), the linear solver should check for solution convergence/accuracy as appropriate.

• SUNLINEARSOLVER_MATRIX_ITERATIVE (2) – the SUNLinSol module requires a matrix, and computes an inexact solution to the linear system defined by that matrix using an iterative algorithm. That is it solves the linear system defined by the matrix object even if that linear system differs from that encoded by the package-supplied ATimes routine. As the solver computes the solution only inexactly (or may diverge), the linear solver should check for solution convergence/accuracy as appropriate.

• SUNLINEARSOLVER_MATRIX_EMBEDDED (3) – the SUNLinSol module sets up and solves the specified linear system at each linear solve call. Any matrix-related data structures are held internally to the linear solver itself, and are not provided by the SUNDIALS package.

Usage:

type = SUNLinSolGetType(LS);

Note

See §11.1.8.1 for more information on intended use cases corresponding to the linear solver type.

SUNLinearSolver_ID SUNLinSolGetID(SUNLinearSolver LS)

Returns a non-negative linear solver identifier (of type int) for the linear solver LS.

Return value:

Non-negative linear solver identifier (of type int), defined by the enumeration SUNLinearSolver_ID, with values shown in Table 11.2 and defined in the sundials_linearsolver.h header file.

Usage:

id = SUNLinSolGetID(LS);

Note

It is recommended that a user-supplied SUNLinearSolver return the SUNLINEARSOLVER_CUSTOM identifier.

SUNErrCode SUNLinSolInitialize(SUNLinearSolver LS)

Performs linear solver initialization (assuming that all solver-specific options have been set).

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolInitialize(LS);

int SUNLinSolSetup(SUNLinearSolver LS, SUNMatrix A)

Performs any linear solver setup needed, based on an updated system SUNMatrix A. This may be called frequently (e.g., with a full Newton method) or infrequently (for a modified Newton method), based on the type of integrator and/or nonlinear solver requesting the solves.

Return value:

Zero for a successful call, a positive value for a recoverable failure, and a negative value for an unrecoverable failure. Ideally this should return one of the generic error codes listed in Table 11.1.

Usage:

retval = SUNLinSolSetup(LS, A);

int SUNLinSolSolve(SUNLinearSolver LS, SUNMatrix A, N_Vector x, N_Vector b, sunrealtype tol)

This required function solves a linear system \(Ax = b\).

Arguments:

• LS – a SUNLinSol object.

• A – a SUNMatrix object.

• x – an N_Vector object containing the initial guess for the solution of the linear system on input, and the solution to the linear system upon return.

• b – an N_Vector object containing the linear system right-hand side.

• tol – the desired linear solver tolerance.

Return value:

Zero for a successful call, a positive value for a recoverable failure, and a negative value for an unrecoverable failure. Ideally this should return one of the generic error codes listed in Table 11.1.

Notes:

Direct solvers: can ignore the tol argument.

Matrix-free solvers: (those that identify as SUNLINEARSOLVER_ITERATIVE) can ignore the SUNMatrix input A, and should rely on the matrix-vector product function supplied through the routine SUNLinSolSetATimes().

Iterative solvers: (those that identify as SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE) should attempt to solve to the specified tolerance tol in a weighted 2-norm. If the solver does not support scaling then it should just use a 2-norm.

Matrix-embedded solvers: should ignore the SUNMatrix input A as this will be NULL. It is assumed that within this function, the solver will call interface routines from the relevant SUNDIALS package to directly form the linear system matrix \(A\), and then solve \(Ax=b\) before returning with the solution \(x\).

Usage:

retval = SUNLinSolSolve(LS, A, x, b, tol);

SUNErrCode SUNLinSolFree(SUNLinearSolver LS)

Frees memory allocated by the linear solver.

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolFree(LS);

11.1.2. SUNLinearSolver “set” functions

The following functions supply linear solver modules with functions defined by the SUNDIALS packages and modify solver parameters. Only the routine for setting the matrix-vector product routine is required, and even then is only required for matrix-free linear solver modules. Otherwise, all other set functions are optional. SUNLinSol implementations that do not provide the functionality for any optional routine should leave the corresponding function pointer NULL instead of supplying a dummy routine.

SUNErrCode SUNLinSolSetATimes(SUNLinearSolver LS, void *A_data, SUNATimesFn ATimes)

Required for matrix-free linear solvers (otherwise optional).

Provides a SUNATimesFn function pointer, as well as a void* pointer to a data structure used by this routine, to the linear solver object LS. SUNDIALS packages call this function to set the matrix-vector product function to either a solver-provided difference-quotient via vector operations or a user-supplied solver-specific routine.

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolSetATimes(LS, A_data, ATimes);

SUNErrCode SUNLinSolSetPreconditioner(SUNLinearSolver LS, void *P_data, SUNPSetupFn Pset, SUNPSolveFn Psol)

This optional routine provides SUNPSetupFn and SUNPSolveFn function pointers that implement the preconditioner solves \(P_1^{-1}\) and \(P_2^{-1}\) from (11.2). This routine is called by a SUNDIALS package, which provides translation between the generic Pset and Psol calls and the package- or user-supplied routines.

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolSetPreconditioner(LS, Pdata, Pset, Psol);

SUNErrCode SUNLinSolSetScalingVectors(SUNLinearSolver LS, N_Vector s1, N_Vector s2)

This optional routine provides left/right scaling vectors for the linear system solve. Here, s1 and s2 are N_Vectors of positive scale factors containing the diagonal of the matrices \(S_1\) and \(S_2\) from (11.2), respectively. Neither vector needs to be tested for positivity, and a NULL argument for either indicates that the corresponding scaling matrix is the identity.

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolSetScalingVectors(LS, s1, s2);

SUNErrCode SUNLinSolSetZeroGuess(SUNLinearSolver LS, sunbooleantype onoff)

This optional routine indicates if the upcoming SUNlinSolSolve() call will be made with a zero initial guess (SUNTRUE) or a non-zero initial guess (SUNFALSE).

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolSetZeroGuess(LS, onoff);

Notes:

It is assumed that the initial guess status is not retained across calls to SUNLinSolSolve(). As such, the linear solver interfaces in each of the SUNDIALS packages call SUNLinSolSetZeroGuess() prior to each call to SUNLinSolSolve().

11.1.3. SUNLinearSolver “get” functions

The following functions allow SUNDIALS packages to retrieve results from a linear solve. All routines are optional.

int SUNLinSolNumIters(SUNLinearSolver LS)

This optional routine should return the number of linear iterations performed in the most-recent “solve” call.

Usage:

its = SUNLinSolNumIters(LS);

sunrealtype SUNLinSolResNorm(SUNLinearSolver LS)

This optional routine should return the final residual norm from the most-recent “solve” call.

Usage:

rnorm = SUNLinSolResNorm(LS);

N_Vector SUNLinSolResid(SUNLinearSolver LS)

If an iterative method computes the preconditioned initial residual and returns with a successful solve without performing any iterations (i.e., either the initial guess or the preconditioner is sufficiently accurate), then this optional routine may be called by the SUNDIALS package. This routine should return the N_Vector containing the preconditioned initial residual vector.

Usage:

rvec = SUNLinSolResid(LS);

Notes:

Since N_Vector is actually a pointer, and the results are not modified, this routine should not require additional memory allocation. If the SUNLinSol object does not retain a vector for this purpose, then this function pointer should be set to NULL in the implementation.

sunindextype SUNLinSolLastFlag(SUNLinearSolver LS)

This optional routine should return the last error flag encountered within the linear solver. Although not called by the SUNDIALS packages directly, this may be called by the user to investigate linear solver issues after a failed solve.

Usage:

lflag = SUNLinLastFlag(LS);

SUNErrCode SUNLinSolSpace(SUNLinearSolver LS, long int *lenrwLS, long int *leniwLS)

This optional routine should return the storage requirements for the linear solver LS:

• lrw is a long int containing the number of sunrealtype words

• liw is a long int containing the number of integer words.

This function is advisory only, for use by users to help determine their total space requirements.

Return value:

A SUNErrCode.

Usage:

retval = SUNLinSolSpace(LS, &lrw, &liw);

11.1.4. Functions provided by SUNDIALS packages

To interface with SUNLinSol modules, the SUNDIALS packages supply a variety of routines for evaluating the matrix-vector product, and setting up and applying the preconditioner. These package-provided routines translate between the user-supplied ODE, DAE, or nonlinear systems and the generic linear solver API. The function types for these routines are defined in the header file sundials/sundials_iterative.h, and are described below.

typedef int (*SUNATimesFn)(void *A_data, N_Vector v, N_Vector z)

Computes the action of a matrix on a vector, performing the operation \(z \gets Av\). Memory for z will already be allocated prior to calling this function. The parameter A_data is a pointer to any information about \(A\) which the function needs in order to do its job. The vector \(v\) should be left unchanged.

Return value:

Zero for a successful call, and non-zero upon failure.

typedef int (*SUNPSetupFn)(void *P_data)

Sets up any requisite problem data in preparation for calls to the corresponding SUNPSolveFn.

Return value:

Zero for a successful call, and non-zero upon failure.

typedef int (*SUNPSolveFn)(void *P_data, N_Vector r, N_Vector z, sunrealtype tol, int lr)

Solves the preconditioner equation \(Pz = r\) for the vector \(z\). Memory for z will already be allocated prior to calling this function. The parameter P_data is a pointer to any information about \(P\) which the function needs in order to do its job (set up by the corresponding SUNPSetupFn). The parameter lr is input, and indicates whether \(P\) is to be taken as the left or right preconditioner: lr = 1 for left and lr = 2 for right. If preconditioning is on one side only, lr can be ignored. If the preconditioner is iterative, then it should strive to solve the preconditioner equation so that

\[\| Pz - r \|_{\text{wrms}} < tol\]

where the error weight vector for the WRMS norm may be accessed from the main package memory structure. The vector r should not be modified by the SUNPSolveFn.

Return value:

Zero for a successful call, a negative value for an unrecoverable failure condition, or a positive value for a recoverable failure condition (thus the calling routine may reattempt the solution after updating preconditioner data).

11.1.5. SUNLinearSolver return codes

The functions provided to SUNLinSol modules by each SUNDIALS package, and functions within the SUNDIALS-provided SUNLinSol implementations, utilize a common set of return codes, listed in Table 11.1. These adhere to a common pattern:

• 0 indicates success

• a positive value corresponds to a recoverable failure, and

• a negative value indicates a non-recoverable failure.

Aside from this pattern, the actual values of each error code provide additional information to the user in case of a linear solver failure.

Table 11.1 SUNLinSol error codes

Error code	Value	Meaning
SUN_SUCCESS	0	successful call or converged solve
SUNLS_ATIMES_NULL	-804	the Atimes function is NULL
SUNLS_ATIMES_FAIL_UNREC	-805	an unrecoverable failure occurred in the ATimes routine
SUNLS_PSET_FAIL_UNREC	-806	an unrecoverable failure occurred in the Pset routine
SUNLS_PSOLVE_NULL	-807	the preconditioner solve function is NULL
SUNLS_PSOLVE_FAIL_UNREC	-808	an unrecoverable failure occurred in the Psolve routine
SUNLS_GS_FAIL	-810	a failure occurred during Gram-Schmidt orthogonalization (SPGMR/SPFGMR)
SUNLS_QRSOL_FAIL	-811	a singular $R$ matrix was encountered in a QR factorization (SPGMR/SPFGMR)
SUNLS_RES_REDUCED	801	an iterative solver reduced the residual, but did not converge to the desired tolerance
SUNLS_CONV_FAIL	802	an iterative solver did not converge (and the residual was not reduced)
SUNLS_ATIMES_FAIL_REC	803	a recoverable failure occurred in the ATimes routine
SUNLS_PSET_FAIL_REC	804	a recoverable failure occurred in the Pset routine
SUNLS_PSOLVE_FAIL_REC	805	a recoverable failure occurred in the Psolve routine
SUNLS_PACKAGE_FAIL_REC	806	a recoverable failure occurred in an external linear solver package
SUNLS_QRFACT_FAIL	807	a singular matrix was encountered during a QR factorization (SPGMR/SPFGMR)
SUNLS_LUFACT_FAIL	808	a singular matrix was encountered during a LU factorization

11.1.6. The generic SUNLinearSolver module

SUNDIALS packages interact with specific SUNLinSol implementations through the generic SUNLinearSolver abstract base class. The SUNLinearSolver type is a pointer to a structure containing an implementation-dependent content field, and an ops field, and is defined as

typedef struct _generic_SUNLinearSolver *SUNLinearSolver

and the generic structure is defined as

struct _generic_SUNLinearSolver {
  void *content;
  struct _generic_SUNLinearSolver_Ops *ops;
};

where the _generic_SUNLinearSolver_Ops structure is a list of pointers to the various actual linear solver operations provided by a specific implementation. The _generic_SUNLinearSolver_Ops structure is defined as

struct _generic_SUNLinearSolver_Ops {
  SUNLinearSolver_Type (*gettype)(SUNLinearSolver);
  SUNLinearSolver_ID   (*getid)(SUNLinearSolver);
  SUNErrCode           (*setatimes)(SUNLinearSolver, void*, SUNATimesFn);
  SUNErrCode           (*setpreconditioner)(SUNLinearSolver, void*,
                                            SUNPSetupFn, SUNPSolveFn);
  SUNErrCode           (*setscalingvectors)(SUNLinearSolver,
                                            N_Vector, N_Vector);
  SUNErrCode           (*setzeroguess)(SUNLinearSolver, sunbooleantype);
  SUNErrCode           (*initialize)(SUNLinearSolver);
  int                  (*setup)(SUNLinearSolver, SUNMatrix);
  int                  (*solve)(SUNLinearSolver, SUNMatrix, N_Vector,
                                N_Vector, sunrealtype);
  int                  (*numiters)(SUNLinearSolver);
  sunrealtype          (*resnorm)(SUNLinearSolver);
  sunindextype         (*lastflag)(SUNLinearSolver);
  SUNErrCode           (*space)(SUNLinearSolver, long int*, long int*);
  N_Vector             (*resid)(SUNLinearSolver);
  SUNErrCode           (*free)(SUNLinearSolver);
};

The generic SUNLinSol class defines and implements the linear solver operations defined in §11.1.1 – §11.1.3. These routines are in fact only wrappers to the linear solver operations defined by a particular SUNLinSol implementation, which are accessed through the ops field of the SUNLinearSolver structure. To illustrate this point we show below the implementation of a typical linear solver operation from the SUNLinearSolver base class, namely SUNLinSolInitialize(), that initializes a SUNLinearSolver object for use after it has been created and configured, and returns a flag denoting a successful or failed operation:

int SUNLinSolInitialize(SUNLinearSolver S)
{
  return ((int) S->ops->initialize(S));
}

11.1.7. Compatibility of SUNLinearSolver modules

Not all SUNLinearSolver implementations are compatible with all SUNMatrix and N_Vector implementations provided in SUNDIALS. More specifically, all of the SUNDIALS iterative linear solvers (SPGMR, SPFGMR, SPBCGS, SPTFQMR, and PCG) are compatible with all of the SUNDIALS N_Vector modules, but the matrix-based direct SUNLinSol modules are specifically designed to work with distinct SUNMatrix and N_Vector modules. In the list below, we summarize the compatibility of each matrix-based SUNLinearSolver module with the various SUNMatrix and N_Vector modules. For a more thorough discussion of these compatibilities, we defer to the documentation for each individual SUNLinSol module in the sections that follow.

• Dense

  • SUNMatrix: Dense or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• LapackDense

  • SUNMatrix: Dense or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• Band

  • SUNMatrix: Band or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• LapackBand

  • SUNMatrix: Band or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• KLU

  • SUNMatrix: Sparse or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• SuperLU_MT

  • SUNMatrix: Sparse or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, or user-supplied

• SuperLU_Dist

  • SUNMatrix: SLUNRLOC or user-supplied

  • N_Vector: Serial, OpenMP, Pthreads, Parallel, *hypre*, PETSc, or user-supplied

• Magma Dense

  • SUNMatrix: Magma Dense or user-supplied

  • N_Vector: HIP, RAJA, or user-supplied

• OneMKL Dense

  • SUNMatrix: One MKL Dense or user-supplied

  • N_Vector: SYCL, RAJA, or user-supplied

• cuSolverSp batchQR

  • SUNMatrix: cuSparse or user-supplied

  • N_Vector: CUDA, RAJA, or user-supplied

11.1.8. Implementing a custom SUNLinearSolver module

A particular implementation of the SUNLinearSolver module must:

• Specify the content field of the SUNLinSol module.

• Define and implement the required linear solver operations.

  Note

  The names of these routines should be unique to that implementation in order to permit using more than one SUNLinSol module (each with different SUNLinearSolver internal data representations) in the same code.

• Define and implement user-callable constructor and destructor routines to create and free a SUNLinearSolver with the new content field and with ops pointing to the new linear solver operations.

We note that the function pointers for all unsupported optional routines should be set to NULL in the ops structure. This allows the SUNDIALS package that is using the SUNLinSol object to know whether the associated functionality is supported.

To aid in the creation of custom SUNLinearSolver modules the generic SUNLinearSolver module provides the utility function SUNLinSolNewEmpty(). When used in custom SUNLinearSolver constructors this function will ease the introduction of any new optional linear solver operations to the SUNLinearSolver API by ensuring that only required operations need to be set.

SUNLinearSolver SUNLinSolNewEmpty()

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

Return value:

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

void SUNLinSolFreeEmpty(SUNLinearSolver LS)

This routine frees the generic SUNLinearSolver 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:

• LS – a SUNLinearSolver object

Additionally, a SUNLinearSolver implementation may do the following:

• Define and implement additional user-callable “set” routines acting on the SUNLinearSolver, e.g., for setting various configuration options to tune the linear solver for a particular problem.

• Provide additional user-callable “get” routines acting on the SUNLinearSolver object, e.g., for returning various solve statistics.

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

Table 11.2 Identifiers associated with SUNLinearSolver modules supplied with SUNDIALS

SUNLinSol ID	Linear solver type	ID Value
SUNLINEARSOLVER_BAND	Banded direct linear solver (internal)	0
SUNLINEARSOLVER_DENSE	Dense direct linear solver (internal)	1
SUNLINEARSOLVER_KLU	Sparse direct linear solver (KLU)	2
SUNLINEARSOLVER_LAPACKBAND	Banded direct linear solver (LAPACK)	3
SUNLINEARSOLVER_LAPACKDENSE	Dense direct linear solver (LAPACK)	4
SUNLINEARSOLVER_PCG	Preconditioned conjugate gradient iterative solver	5
SUNLINEARSOLVER_SPBCGS	Scaled-preconditioned BiCGStab iterative solver	6
SUNLINEARSOLVER_SPFGMR	Scaled-preconditioned FGMRES iterative solver	7
SUNLINEARSOLVER_SPGMR	Scaled-preconditioned GMRES iterative solver	8
SUNLINEARSOLVER_SPTFQMR	Scaled-preconditioned TFQMR iterative solver	9
SUNLINEARSOLVER_SUPERLUDIST	Parallel sparse direct linear solver (SuperLU_Dist)	10
SUNLINEARSOLVER_SUPERLUMT	Threaded sparse direct linear solver (SuperLU_MT)	11
SUNLINEARSOLVER_CUSOLVERSP_BATCHQR	Sparse direct linear solver (CUDA)	12
SUNLINEARSOLVER_MAGMADENSE	Dense or block-dense direct linear solver (MAGMA)	13
SUNLINEARSOLVER_ONEMKLDENSE	Dense or block-dense direct linear solver (OneMKL)	14
SUNLINEARSOLVER_CUSTOM	User-provided custom linear solver	15

11.1.8.1. Intended use cases

The SUNLinSol and SUNMATRIX APIs are designed to require a minimal set of routines to ease interfacing with custom or third-party linear solver libraries. Many external solvers provide routines with similar functionality and thus may require minimal effort to wrap within custom SUNMATRIX and SUNLinSol implementations. As SUNDIALS packages utilize generic SUNLinSol modules they may naturally leverage user-supplied SUNLinearSolver implementations, thus there exist a wide range of possible linear solver combinations. Some intended use cases for both the SUNDIALS-provided and user-supplied SUNLinSol modules are discussd in the sections below.

11.1.8.1.1. Direct linear solvers

Direct linear solver modules require a matrix and compute an “exact” solution to the linear system defined by the matrix. SUNDIALS packages strive to amortize the high cost of matrix construction by reusing matrix information for multiple nonlinear iterations or time steps. As a result, each package’s linear solver interface recomputes matrix information as infrequently as possible.

Alternative matrix storage formats and compatible linear solvers that are not currently provided by, or interfaced with, SUNDIALS can leverage this infrastructure with minimal effort. To do so, a user must implement custom SUNMATRIX and SUNLinSol wrappers for the desired matrix format and/or linear solver following the APIs described in §10 and §11. This user-supplied SUNLinSol module must then self-identify as having SUNLINEARSOLVER_DIRECT type.

11.1.8.1.2. Matrix-free iterative linear solvers

Matrix-free iterative linear solver modules do not require a matrix, and instead compute an inexact solution to the linear system defined by the package-supplied ATimes routine. SUNDIALS supplies multiple scaled, preconditioned iterative SUNLinSol modules that support scaling, allowing packages to handle non-dimensionalization, and users to define variables and equations as natural in their applications. However, for linear solvers that do not support left/right scaling, SUNDIALS packages must instead adjust the tolerance supplied to the linear solver to compensate (see the iterative linear tolerance section that follows for more details) – this strategy may be non-optimal since it cannot handle situations where the magnitudes of different solution components or equations vary dramatically within a single application.

To utilize alternative linear solvers that are not currently provided by, or interfaced with, SUNDIALS a user must implement a custom SUNLinSol wrapper for the linear solver following the API described in §11. This user-supplied SUNLinSol module must then self-identify as having SUNLINEARSOLVER_ITERATIVE type.

11.1.8.1.3. Matrix-based iterative linear solvers (reusing \(A\))

Matrix-based iterative linear solver modules require a matrix and compute an inexact solution to the linear system defined by the matrix. This matrix will be updated infrequently and resued across multiple solves to amortize the cost of matrix construction. As in the direct linear solver case, only thin SUNMATRIX and SUNLinSol wrappers for the underlying matrix and linear solver structures need to be created to utilize such a linear solver. This user-supplied SUNLinSol module must then self-identify as having SUNLINEARSOLVER_MATRIX_ITERATIVE type.

At present, SUNDIALS has one example problem that uses this approach for wrapping a structured-grid matrix, linear solver, and preconditioner from the hypre library; this may be used as a template for other customized implementations (see examples/arkode/CXX_parhyp/ark_heat2D_hypre.cpp).

11.1.8.1.4. Matrix-based iterative linear solvers (current \(A\))

For users who wish to utilize a matrix-based iterative linear solver where the matrix is purely for preconditioning and the linear system is defined by the package-supplied ATimes routine, we envision two current possibilities.

The preferred approach is for users to employ one of the SUNDIALS scaled, preconditioned iterative linear solver implementations (SUNLinSol_SPGMR(), SUNLinSol_SPFGMR(), SUNLinSol_SPBCGS(), SUNLinSol_SPTFQMR(), or SUNLinSol_PCG()) as the outer solver. The creation and storage of the preconditioner matrix, and interfacing with the corresponding matrix-based linear solver, can be handled through a package’s preconditioner “setup” and “solve” functionality without creating SUNMATRIX and SUNLinSol implementations. This usage mode is recommended primarily because the SUNDIALS-provided modules support variable and equation scaling as described above.

A second approach supported by the linear solver APIs is as follows. If the SUNLinSol implementation is matrix-based, self-identifies as having SUNLINEARSOLVER_ITERATIVE type, and also provides a non-NULL SUNLinSolSetATimes() routine, then each SUNDIALS package will call that routine to attach its package-specific matrix-vector product routine to the SUNLinSol object. The SUNDIALS package will then call the SUNLinSol-provided SUNLinSolSetup() routine (infrequently) to update matrix information, but will provide current matrix-vector products to the SUNLinSol implementation through the package-supplied SUNATimesFn routine.

11.1.8.1.5. Application-specific linear solvers with embedded matrix structure

Many applications can exploit additional linear system structure arising from to the implicit couplings in their model equations. In certain circumstances, the linear solve \(Ax=b\) may be performed without the need for a global system matrix \(A\), as the unformed \(A\) may be block diagonal or block triangular, and thus the overall linear solve may be performed through a sequence of smaller linear solves. In other circumstances, a linear system solve may be accomplished via specialized fast solvers, such as the fast Fourier transform, fast multipole method, or treecode, in which case no matrix structure may be explicitly necessary. In many of the above situations, construction and preprocessing of the linear system matrix \(A\) may be inexpensive, and thus increased performance may be possible if the current linear system information is used within every solve (instead of being lagged, as occurs with matrix-based solvers that reuse \(A\)).

To support such application-specific situations, SUNDIALS supports user-provided linear solvers with the SUNLINEARSOLVER_MATRIX_EMBEDDED type. For an application to leverage this support, it should define a custom SUNLinSol implementation having this type, that only needs to implement the required SUNLinSolGetType() and SUNLinSolSolve() operations. Within SUNLinSolSolve(), the linear solver implementation should call package-specific interface routines (e.g., ARKStepGetNonlinearSystemData, CVodeGetNonlinearSystemData, IDAGetNonlinearSystemData, ARKStepGetCurrentGamma, CVodeGetCurrentGamma, IDAGetCurrentCj, or MRIStepGetCurrentGamma) to construct the relevant system matrix \(A\) (or portions thereof), solve the linear system \(Ax=b\), and return the solution vector \(x\).

We note that when attaching this custom SUNLinearSolver object with the relevant SUNDIALS package SetLinearSolver routine, the input SUNMatrix A should be set to NULL.

For templates of such user-provided “matrix-embedded” SUNLinSol implementations, see the SUNDIALS examples ark_analytic_mels.c, cvAnalytic_mels.c, cvsAnalytic_mels.c, idaAnalytic_mels.c, and idasAnalytic_mels.c.