10.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
.
10.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 solverspecific 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()
).

enum SUNLinearSolver_Type
An identifier indicating the type of linear solver.
Note
See §10.1.8.1 for more information on intended use cases corresponding to the linear solver type.

enumerator SUNLINEARSOLVER_DIRECT
The linear solver requires a matrix, and computes an “exact” solution to the linear system defined by that matrix.

enumerator SUNLINEARSOLVER_ITERATIVE
The linear solver does not require a matrix (though one may be provided), and computes an inexact solution to the linear system using a matrixfree iterative algorithm. That is it solves the linear system defined by the packagesupplied
ATimes
routine (seeSUNLinSolSetATimes()
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.

enumerator SUNLINEARSOLVER_MATRIX_ITERATIVE
The linear solver 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 packagesupplied
ATimes
routine. As the solver computes the solution only inexactly (or may diverge), the linear solver should check for solution convergence/accuracy as appropriate.

enumerator SUNLINEARSOLVER_MATRIX_EMBEDDED
The linear solver sets up and solves the specified linear system at each linear solve call. Any matrixrelated data structures are held internally to the linear solver itself, and are not provided by the SUNDIALS package.

enumerator SUNLINEARSOLVER_DIRECT

SUNLinearSolver_Type SUNLinSolGetType(SUNLinearSolver LS)
Returns the
SUNLinearSolver_Type
type identifier for the linear solver.Usage:
type = SUNLinSolGetType(LS);

SUNLinearSolver_ID SUNLinSolGetID(SUNLinearSolver LS)
Returns a nonnegative linear solver identifier (of type
int
) for the linear solver LS.Return value:
Nonnegative linear solver identifier (of type
int
), defined by the enumerationSUNLinearSolver_ID
, with values shown in Table 10.2 and defined in thesundials_linearsolver.h
header file.Usage:
id = SUNLinSolGetID(LS);
Note
It is recommended that a usersupplied
SUNLinearSolver
return theSUNLINEARSOLVER_CUSTOM
identifier.

SUNErrCode SUNLinSolInitialize(SUNLinearSolver LS)
Performs linear solver initialization (assuming that all solverspecific 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 10.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 righthand 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 10.1.
Notes:
Direct solvers: can ignore the tol argument.
Matrixfree solvers: (those that identify as
SUNLINEARSOLVER_ITERATIVE
) can ignore theSUNMatrix
input A, and should rely on the matrixvector product function supplied through the routineSUNLinSolSetATimes()
.Iterative solvers: (those that identify as
SUNLINEARSOLVER_ITERATIVE
orSUNLINEARSOLVER_MATRIX_ITERATIVE
) should attempt to solve to the specified tolerance tol in a weighted 2norm. If the solver does not support scaling then it should just use a 2norm.Matrixembedded solvers: should ignore the
SUNMatrix
input A as this will beNULL
. 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);
10.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 matrixvector product routine is required, and even then is
only required for matrixfree 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 matrixfree linear solvers (otherwise optional).
Provides a
SUNATimesFn
function pointer, as well as avoid*
pointer to a data structure used by this routine, to the linear solver object LS. SUNDIALS packages call this function to set the matrixvector product function to either a solverprovided differencequotient via vector operations or a usersupplied solverspecific 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
andSUNPSolveFn
function pointers that implement the preconditioner solves \(P_1^{1}\) and \(P_2^{1}\) from (10.2). This routine is called by a SUNDIALS package, which provides translation between the generic Pset and Psol calls and the package or usersupplied 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 vectors of positive scale factors containing the diagonal of the matrices \(S_1\) and \(S_2\) from (10.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 nonzero 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 callSUNLinSolSetZeroGuess()
prior to each call toSUNLinSolSolve()
.
10.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 mostrecent “solve” call.
Usage:
its = SUNLinSolNumIters(LS);

sunrealtype SUNLinSolResNorm(SUNLinearSolver LS)
This optional routine should return the final residual norm from the mostrecent “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 toNULL
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 wordsliw 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);
10.1.4. Functions provided by SUNDIALS packages
To interface with SUNLinSol modules, the SUNDIALS packages supply a
variety of routines for evaluating the matrixvector product, and
setting up and applying the preconditioner. These packageprovided
routines translate between the usersupplied 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 nonzero 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 nonzero 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).
10.1.5. SUNLinearSolver return codes
The functions provided to SUNLinSol modules by each SUNDIALS package, and functions within the SUNDIALSprovided SUNLinSol implementations, utilize a common set of return codes, listed in Table 10.1. These adhere to a common pattern:
0 indicates success
a positive value corresponds to a recoverable failure, and
a negative value indicates a nonrecoverable 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.
Error code 
Value 
Meaning 


0 
successful call or converged solve 

804 
the 

805 
an unrecoverable failure occurred in the


806 
an unrecoverable failure occurred in the 

807 
the preconditioner solve function is 

808 
an unrecoverable failure occurred in the


810 
a failure occurred during GramSchmidt orthogonalization (SPGMR/SPFGMR) 

811 
a singular $R$ matrix was encountered in a QR factorization (SPGMR/SPFGMR) 

801 
an iterative solver reduced the residual, but did not converge to the desired tolerance 

802 
an iterative solver did not converge (and the residual was not reduced) 

803 
a recoverable failure occurred in the 

804 
a recoverable failure occurred in the 

805 
a recoverable failure occurred in the 

806 
a recoverable failure occurred in an external linear solver package 

807 
a singular matrix was encountered during a QR factorization (SPGMR/SPFGMR) 

808 
a singular matrix was encountered during a LU factorization 
10.1.6. The generic SUNLinearSolver module
SUNDIALS packages interact with linear solver implementations through the
SUNLinearSolver
class. A SUNLinearSolver
is a pointer to the
_generic_SUNLinearSolver
structure:

typedef struct _generic_SUNLinearSolver *SUNLinearSolver

struct _generic_SUNLinearSolver
The structure defining the SUNDIALS linear solver class.

void *content
Pointer to the linear solverspecific member data

SUNLinearSolver_Ops ops
A virtual table of linear solver operations provided by a specific implementation

SUNContext sunctx
The SUNDIALS simulation context

void *content
The virtual table structure is defined as

typedef struct _generic_SUNLinearSolver_Ops *SUNLinearSolver_Ops

struct _generic_SUNLinearSolver_Ops
The structure defining
SUNLinearSolver
operations.
SUNLinearSolver_Type (*gettype)(SUNLinearSolver)
The function implementing
SUNLinSolGetType()

SUNLinearSolver_ID (*getid)(SUNLinearSolver)
The function implementing
SUNLinSolGetID()

SUNErrCode (*setatimes)(SUNLinearSolver, void*, SUNATimesFn)
The function implementing
SUNLinSolSetATimes()

SUNErrCode (*setpreconditioner)(SUNLinearSolver, void*, SUNPSetupFn, SUNPSolveFn)
The function implementing
SUNLinSolSetPreconditioner()

SUNErrCode (*setscalingvectors)(SUNLinearSolver, N_Vector, N_Vector)
The function implementing
SUNLinSolSetScalingVectors()

SUNErrCode (*setzeroguess)(SUNLinearSolver, sunbooleantype)
The function implementing
SUNLinSolSetZeroGuess()

SUNErrCode (*initialize)(SUNLinearSolver)
The function implementing
SUNLinSolInitialize()

int (*setup)(SUNLinearSolver, SUNMatrix)
The function implementing
SUNLinSolSetup()

int (*solve)(SUNLinearSolver, SUNMatrix, N_Vector, N_Vector, sunrealtype)
The function implementing
SUNLinSolSolve()

int (*numiters)(SUNLinearSolver)
The function implementing
SUNLinSolNumIters()

sunrealtype (*resnorm)(SUNLinearSolver)
The function implementing
SUNLinSolResNorm()

sunindextype (*lastflag)(SUNLinearSolver)
The function implementing
SUNLinSolLastFlag()

SUNErrCode (*space)(SUNLinearSolver, long int*, long int*)
The function implementing
SUNLinSolSpace()

N_Vector (*resid)(SUNLinearSolver)
The function implementing
SUNLinSolResid()

SUNErrCode (*free)(SUNLinearSolver)
The function implementing
SUNLinSolFree()

SUNLinearSolver_Type (*gettype)(SUNLinearSolver)
The generic SUNLinSol class defines and implements the linear solver
operations defined in §10.1.1 – §10.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));
}
10.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 matrixbased 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 matrixbased 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.

SUNMatrix
: Magma Dense or usersupplied

SUNMatrix
: One MKL Dense or usersupplied
10.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 usercallable 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(SUNContext sunctx)
This function allocates a new generic
SUNLinearSolver
object and initializes its content pointer and the function pointers in the operations structure toNULL
.Return value:
If successful, this function returns a
SUNLinearSolver
object. If an error occurs when allocating the object, then this routine will returnNULL
.

void SUNLinSolFreeEmpty(SUNLinearSolver LS)
This routine frees the generic
SUNLinearSolver
object, under the assumption that any implementationspecific data that was allocated within the underlying content structure has already been freed. It will additionally test whether the ops pointer isNULL
, 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 usercallable “set” routines acting on the
SUNLinearSolver
, e.g., for setting various configuration options to tune the linear solver for a particular problem.Provide additional usercallable “get” routines acting on the
SUNLinearSolver
object, e.g., for returning various solve statistics.

enum SUNLinearSolver_ID
Each SUNLinSol implementation included in SUNDIALS has a unique identifier specified in enumeration and shown in Table 10.2. It is recommended that a usersupplied SUNLinSol implementation use the
SUNLINEARSOLVER_CUSTOM
identifier.
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 
Scaledpreconditioned BiCGStab iterative solver 
6 
SUNLINEARSOLVER_SPFGMR 
Scaledpreconditioned FGMRES iterative solver 
7 
SUNLINEARSOLVER_SPGMR 
Scaledpreconditioned GMRES iterative solver 
8 
SUNLINEARSOLVER_SPTFQMR 
Scaledpreconditioned 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 blockdense direct linear solver (MAGMA) 
13 
SUNLINEARSOLVER_ONEMKLDENSE 
Dense or blockdense direct linear solver (OneMKL) 
14 
SUNLINEARSOLVER_CUSTOM 
Userprovided custom linear solver 
15 
10.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 thirdparty 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 usersupplied
SUNLinearSolver
implementations, thus there exist a wide range of
possible linear solver combinations. Some intended use cases for both the
SUNDIALSprovided and usersupplied SUNLinSol modules are discussd in the
sections below.
10.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 §9
and §10. This usersupplied SUNLinSol module must then
selfidentify as having SUNLINEARSOLVER_DIRECT
type.
10.1.8.1.2. Matrixfree iterative linear solvers
Matrixfree iterative linear solver modules do not require a matrix, and instead
compute an inexact solution to the linear system defined by the
packagesupplied ATimes
routine. SUNDIALS supplies multiple scaled,
preconditioned iterative SUNLinSol modules that support scaling, allowing
packages to handle nondimensionalization, 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
nonoptimal 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
§10. This usersupplied SUNLinSol module must then
selfidentify as having SUNLINEARSOLVER_ITERATIVE
type.
10.1.8.1.3. Matrixbased iterative linear solvers (reusing \(A\))
Matrixbased 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 usersupplied SUNLinSol module must then
selfidentify as having SUNLINEARSOLVER_MATRIX_ITERATIVE
type.
At present, SUNDIALS has one example problem that uses this approach for
wrapping a structuredgrid 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
).
10.1.8.1.4. Matrixbased iterative linear solvers (current \(A\))
For users who wish to utilize a matrixbased iterative linear solver
where the matrix is purely for preconditioning and the linear system is
defined by the packagesupplied 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 matrixbased
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 SUNDIALSprovided 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 matrixbased, selfidentifies
as having SUNLINEARSOLVER_ITERATIVE
type, and also provides a nonNULL
SUNLinSolSetATimes()
routine, then each SUNDIALS package
will call that routine to attach its packagespecific matrixvector
product routine to the SUNLinSol object. The SUNDIALS package will
then call the SUNLinSolprovided SUNLinSolSetup()
routine
(infrequently) to update matrix information, but will provide current
matrixvector products to the SUNLinSol implementation through the
packagesupplied SUNATimesFn
routine.
10.1.8.1.5. Applicationspecific 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 matrixbased solvers that reuse \(A\)).
To support such applicationspecific situations, SUNDIALS supports userprovided
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 packagespecific 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 userprovided “matrixembedded” SUNLinSol implementations,
see the SUNDIALS examples ark_analytic_mels.c
, cvAnalytic_mels.c
,
cvsAnalytic_mels.c
, idaAnalytic_mels.c
, and idasAnalytic_mels.c
.