10.8. The SUNLinSol_Band Module
The SUNLinSol_Band implementation of the SUNLinearSolver
class
is designed to be used with the corresponding
SUNMATRIX_BAND matrix type, and one of the
serial or sharedmemory N_Vector
implementations
(NVECTOR_SERIAL, NVECTOR_OPENMP or NVECTOR_PTHREADS).
10.8.1. SUNLinSol_Band Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_band.h
. The SUNLinSol_Band module
is accessible from all SUNDIALS packages without
linking to the
libsundials_sunlinsolband
module library.
The SUNLinSol_Band module provides the following usercallable constructor routine:

SUNLinearSolver SUNLinSol_Band(N_Vector y, SUNMatrix A, SUNContext sunctx)
This function creates and allocates memory for a band
SUNLinearSolver
. Arguments:
y – vector used to determine the linear system size
A – matrix used to assess compatibility
sunctx – the
SUNContext
object (see §1.4)
 Return value:
New SUNLinSol_Band object, or
NULL
if eitherA
ory
are incompatible. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_BAND matrix type and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.Additionally, this routine will verify that the input matrix
A
is allocated with appropriate upper bandwidth storage for the \(LU\) factorization.
10.8.2. SUNLinSol_Band Description
The SUNLinSol_Band module defines the content
field of a SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_Band {
sunindextype N;
sunindextype *pivots;
sunindextype last_flag;
};
These entries of the content field contain the following information:
N
 size of the linear system,pivots
 index array for partial pivoting in LU factorization,last_flag
 last error return flag from internal function evaluations.
This solver is constructed to perform the following operations:
The “setup” call performs an \(LU\) factorization with partial (row) pivoting, \(PA=LU\), where \(P\) is a permutation matrix, \(L\) is a lower triangular matrix with 1’s on the diagonal, and \(U\) is an upper triangular matrix. This factorization is stored inplace on the input SUNMATRIX_BAND object \(A\), with pivoting information encoding \(P\) stored in the
pivots
array.The “solve” call performs pivoting and forward and backward substitution using the stored
pivots
array and the \(LU\) factors held in the SUNMATRIX_BAND object.\(A\) must be allocated to accommodate the increase in upper bandwidth that occurs during factorization. More precisely, if \(A\) is a band matrix with upper bandwidth
mu
and lower bandwidthml
, then the upper triangular factor \(U\) can have upper bandwidth as big assmu = MIN(N1,mu+ml)
. The lower triangular factor \(L\) has lower bandwidthml
.
The SUNLinSol_Band module defines band implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_Band
SUNLinSolInitialize_Band
– this does nothing, since all consistency checks are performed at solver creation.SUNLinSolSetup_Band
– this performs the \(LU\) factorization.SUNLinSolSolve_Band
– this uses the \(LU\) factors andpivots
array to perform the solve.SUNLinSolLastFlag_Band
SUNLinSolSpace_Band
– this only returns information for the storage within the solver object, i.e. storage forN
,last_flag
, andpivots
.SUNLinSolFree_Band
10.9. The SUNLinSol_Dense Module
The SUNLinSol_Dense implementation of the SUNLinearSolver
class
is designed to be used with the corresponding SUNMATRIX_DENSE matrix type,
and one of the serial or sharedmemory N_Vector
implementations
(NVECTOR_SERIAL, NVECTOR_OPENMP or NVECTOR_PTHREADS).
10.9.1. SUNLinSol_Dense Usage
The header file to be included when using this module is
sunlinsol/sunlinsol_dense.h
. The SUNLinSol_Dense module is
accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsoldense
module library.
The module SUNLinSol_Dense provides the following usercallable constructor routine:

SUNLinearSolver SUNLinSol_Dense(N_Vector y, SUNMatrix A, SUNContext sunctx)
This function creates and allocates memory for a dense
SUNLinearSolver
. Arguments:
y – vector used to determine the linear system size.
A – matrix used to assess compatibility.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
New SUNLinSol_Dense object, or
NULL
if eitherA
ory
are incompatible. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_DENSE matrix type and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.
10.9.2. SUNLinSol_Dense Description
The SUNLinSol_Dense module defines the content
field of a SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_Dense {
sunindextype N;
sunindextype *pivots;
sunindextype last_flag;
};
These entries of the content field contain the following information:
N
 size of the linear system,pivots
 index array for partial pivoting in LU factorization,last_flag
 last error return flag from internal function evaluations.
This solver is constructed to perform the following operations:
The “setup” call performs an \(LU\) factorization with partial (row) pivoting (\(\mathcal O(N^3)\) cost), \(PA=LU\), where \(P\) is a permutation matrix, \(L\) is a lower triangular matrix with 1’s on the diagonal, and \(U\) is an upper triangular matrix. This factorization is stored inplace on the input SUNMATRIX_DENSE object \(A\), with pivoting information encoding \(P\) stored in the
pivots
array.The “solve” call performs pivoting and forward and backward substitution using the stored
pivots
array and the \(LU\) factors held in the SUNMATRIX_DENSE object (\(\mathcal O(N^2)\) cost).
The SUNLinSol_Dense module defines dense implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_Dense
SUNLinSolInitialize_Dense
– this does nothing, since all consistency checks are performed at solver creation.SUNLinSolSetup_Dense
– this performs the \(LU\) factorization.SUNLinSolSolve_Dense
– this uses the \(LU\) factors andpivots
array to perform the solve.SUNLinSolLastFlag_Dense
SUNLinSolSpace_Dense
– this only returns information for the storage within the solver object, i.e. storage forN
,last_flag
, andpivots
.SUNLinSolFree_Dense
10.10. The SUNLinSol_KLU Module
The SUNLinSol_KLU implementation of the SUNLinearSolver
class
is designed to be used with the corresponding SUNMATRIX_SPARSE matrix type,
and one of the serial or sharedmemory N_Vector
implementations
(NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS).
10.10.1. SUNLinSol_KLU Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_klu.h
. The installed module
library to link to is libsundials_sunlinsolklu
.lib
where .lib is typically .so
for shared libraries and
.a
for static libraries.
The module SUNLinSol_KLU provides the following additional usercallable routines:

SUNLinearSolver SUNLinSol_KLU(N_Vector y, SUNMatrix A, SUNContext sunctx)
This constructor function creates and allocates memory for a SUNLinSol_KLU object.
 Arguments:
y – vector used to determine the linear system size.
A – matrix used to assess compatibility.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
New SUNLinSol_KLU object, or
NULL
if eitherA
ory
are incompatible. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_SPARSE matrix type (using either CSR or CSC storage formats) and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.

SUNErrCode SUNLinSol_KLUReInit(SUNLinearSolver S, SUNMatrix A, sunindextype nnz, int reinit_type)
This function reinitializes memory and flags for a new factorization (symbolic and numeric) to be conducted at the next solver setup call. This routine is useful in the cases where the number of nonzeroes has changed or if the structure of the linear system has changed which would require a new symbolic (and numeric factorization).
 Arguments:
S – existing SUNLinSol_KLU object to reinitialize.
A – sparse
SUNMatrix
matrix (with updated structure) to use for reinitialization.nnz – maximum number of nonzeros expected for Jacobian matrix.
reinit_type – governs the level of reinitialization. The allowed values are:
The Jacobian matrix will be destroyed and a new one will be allocated based on the
nnz
value passed to this call. New symbolic and numeric factorizations will be completed at the next solver setup.Only symbolic and numeric factorizations will be completed. It is assumed that the Jacobian size has not exceeded the size of
nnz
given in the sparse matrix provided to the original constructor routine (or the previousSUNKLUReInit
call).
 Return value:
 Notes:
This routine assumes no other changes to solver use are necessary.

SUNErrCode SUNLinSol_KLUSetOrdering(SUNLinearSolver S, int ordering_choice)
This function sets the ordering used by KLU for reducing fill in the linear solve.
 Arguments:
S – existing SUNLinSol_KLU object to update.
ordering_choice – type of ordering to use, options are:
AMD,
COLAMD, and
the natural ordering.
The default is 1 for COLAMD.
 Return value:

sun_klu_symbolic *SUNLinSol_KLUGetSymbolic(SUNLinearSolver S)
This function returns a pointer to the KLU symbolic factorization stored in the SUNLinSol_KLU
content
structure.
type sun_klu_symbolic
This type is an alias that depends on the SUNDIALS index size, see
sunindextype
andSUNDIALS_INDEX_SIZE
. It is equivalent to:klu_symbolic
when SUNDIALS is compiled with 32bit indicesklu_l_symbolic
when SUNDIALS is compiled with 64bit indices

type sun_klu_symbolic

sun_klu_numeric *SUNLinSol_KLUGetNumeric(SUNLinearSolver S)
This function returns a pointer to the KLU numeric factorization stored in the SUNLinSol_KLU
content
structure.
type sun_klu_numeric
This type is an alias that depends on the SUNDIALS index size, see
sunindextype
andSUNDIALS_INDEX_SIZE
. It is equivalent to:klu_numeric
when SUNDIALS is compiled with 32bit indicesklu_l_numeric
when SUNDIALS is compiled with 64bit indices

type sun_klu_numeric

sun_klu_common *SUNLinSol_KLUGetCommon(SUNLinearSolver S)
This function returns a pointer to the KLU common structure stored in the SUNLinSol_KLU
content
structure.
type sun_klu_common
This type is an alias that depends on the SUNDIALS index size, see
sunindextype
andSUNDIALS_INDEX_SIZE
. It is equivalent to:klu_common
when SUNDIALS is compiled with 32bit indicesklu_l_common
when SUNDIALS is compiled with 64bit indices

type sun_klu_common
10.10.2. SUNLinSol_KLU Description
The SUNLinSol_KLU module defines the content
field of a SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_KLU {
int last_flag;
int first_factorize;
sun_klu_symbolic *symbolic;
sun_klu_numeric *numeric;
sun_klu_common common;
sunindextype (*klu_solver)(sun_klu_symbolic*, sun_klu_numeric*,
sunindextype, sunindextype,
double*, sun_klu_common*);
};
These entries of the content field contain the following information:
last_flag
 last error return flag from internal function evaluations,first_factorize
 flag indicating whether the factorization has ever been performed,symbolic
 KLU storage structure for symbolic factorization components, with underlying typeklu_symbolic
orklu_l_symbolic
, depending on whether SUNDIALS was installed with 32bit versus 64bit indices, respectively,numeric
 KLU storage structure for numeric factorization components, with underlying typeklu_numeric
orklu_l_numeric
, depending on whether SUNDIALS was installed with 32bit versus 64bit indices, respectively,common
 storage structure for common KLU solver components, with underlying typeklu_common
orklu_l_common
, depending on whether SUNDIALS was installed with 32bit versus 64bit indices, respectively,klu_solver
– pointer to the appropriate KLU solver function (depending on whether it is using a CSR or CSC sparse matrix, and on whether SUNDIALS was installed with 32bit or 64bit indices).
The SUNLinSol_KLU module is a SUNLinearSolver
wrapper for
the KLU sparse matrix factorization and solver library written by Tim
Davis and collaborators ([4, 40]). In order to use the
SUNLinSol_KLU interface to KLU, it is assumed that KLU has
been installed on the system prior to installation of SUNDIALS, and
that SUNDIALS has been configured appropriately to link with KLU
(see §1.2.4 for details).
Additionally, this wrapper only supports doubleprecision
calculations, and therefore cannot be compiled if SUNDIALS is
configured to have sunrealtype
set to either extended
or
single
(see §1.3 for
details). Since the KLU library supports both 32bit and 64bit
integers, this interface will be compiled for either of the available
sunindextype
options.
The KLU library has a symbolic factorization routine that computes the permutation of the linear system matrix to block triangular form and the permutations that will preorder the diagonal blocks (the only ones that need to be factored) to reduce fillin (using AMD, COLAMD, CHOLAMD, natural, or an ordering given by the user). Of these ordering choices, the default value in the SUNLinSol_KLU module is the COLAMD ordering.
KLU breaks the factorization into two separate parts. The first is a symbolic factorization and the second is a numeric factorization that returns the factored matrix along with final pivot information. KLU also has a refactor routine that can be called instead of the numeric factorization. This routine will reuse the pivot information. This routine also returns diagnostic information that a user can examine to determine if numerical stability is being lost and a full numerical factorization should be done instead of the refactor.
Since the linear systems that arise within the context of SUNDIALS calculations will typically have identical sparsity patterns, the SUNLinSol_KLU module is constructed to perform the following operations:
The first time that the “setup” routine is called, it performs the symbolic factorization, followed by an initial numerical factorization.
On subsequent calls to the “setup” routine, it calls the appropriate KLU “refactor” routine, followed by estimates of the numerical conditioning using the relevant “rcond”, and if necessary “condest”, routine(s). If these estimates of the condition number are larger than \(\varepsilon^{2/3}\) (where \(\varepsilon\) is the doubleprecision unit roundoff), then a new factorization is performed.
The module includes the routine
SUNKLUReInit
, that can be called by the user to force a full refactorization at the next “setup” call.The “solve” call performs pivoting and forward and backward substitution using the stored KLU data structures. We note that in this solve KLU operates on the native data arrays for the righthand side and solution vectors, without requiring costly data copies.
The SUNLinSol_KLU module defines implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_KLU
SUNLinSolInitialize_KLU
– this sets thefirst_factorize
flag to 1, forcing both symbolic and numerical factorizations on the subsequent “setup” call.SUNLinSolSetup_KLU
– this performs either a \(LU\) factorization or refactorization of the input matrix.SUNLinSolSolve_KLU
– this calls the appropriate KLU solve routine to utilize the \(LU\) factors to solve the linear system.SUNLinSolLastFlag_KLU
SUNLinSolSpace_KLU
– this only returns information for the storage within the solver interface, i.e. storage for the integerslast_flag
andfirst_factorize
. For additional space requirements, see the KLU documentation.SUNLinSolFree_KLU
10.11. The SUNLinSol_LapackBand Module
The SUNLinSol_LapackBand implementation of the SUNLinearSolver
class
is designed to be used with the corresponding SUNMATRIX_BAND matrix type,
and one of the serial or sharedmemory N_Vector
implementations
(NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS). The
10.11.1. SUNLinSol_LapackBand Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_lapackband.h
. The installed module
library to link to is
libsundials_sunlinsollapackband
.lib
where .lib is typically .so
for shared libraries and
.a
for static libraries.
The module SUNLinSol_LapackBand provides the following usercallable routine:

SUNLinearSolver SUNLinSol_LapackBand(N_Vector y, SUNMatrix A, SUNContext sunctx)
This function creates and allocates memory for a LAPACK band
SUNLinearSolver
. Arguments:
y – vector used to determine the linear system size.
A – matrix used to assess compatibility.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
New SUNLinSol_LapackBand object, or
NULL
if eitherA
ory
are incompatible. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_BAND matrix type and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.Additionally, this routine will verify that the input matrix
A
is allocated with appropriate upper bandwidth storage for the \(LU\) factorization.
10.11.2. SUNLinSol_LapackBand Description
SUNLinSol_LapackBand module defines the
content field of a SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_Band {
sunindextype N;
sunindextype *pivots;
sunindextype last_flag;
};
These entries of the content field contain the following information:
N
 size of the linear system,pivots
 index array for partial pivoting in LU factorization,last_flag
 last error return flag from internal function evaluations.
The SUNLinSol_LapackBand module is a SUNLinearSolver
wrapper for
the LAPACK band matrix factorization and solve routines, *GBTRF
and *GBTRS
, where *
is either D
or S
, depending on
whether SUNDIALS was configured to have sunrealtype
set to
double
or single
, respectively (see
§1.3 for details).
In order to use the SUNLinSol_LapackBand module it is assumed
that LAPACK has been installed on the system prior to installation of
SUNDIALS, and that SUNDIALS has been configured appropriately to
link with LAPACK (see
§1.2.4 for details). We note
that since there do not exist 128bit floatingpoint factorization and
solve routines in LAPACK, this interface cannot be compiled when
using extended
precision for sunrealtype
. Similarly, since
there do not exist 64bit integer LAPACK routines, the
SUNLinSol_LapackBand module also cannot be compiled when using
int64_t
for the sunindextype
.
This solver is constructed to perform the following operations:
The “setup” call performs an \(LU\) factorization with partial (row) pivoting, \(PA=LU\), where \(P\) is a permutation matrix, \(L\) is a lower triangular matrix with 1’s on the diagonal, and \(U\) is an upper triangular matrix. This factorization is stored inplace on the input SUNMATRIX_BAND object \(A\), with pivoting information encoding \(P\) stored in the
pivots
array.The “solve” call performs pivoting and forward and backward substitution using the stored
pivots
array and the \(LU\) factors held in the SUNMATRIX_BAND object.\(A\) must be allocated to accommodate the increase in upper bandwidth that occurs during factorization. More precisely, if \(A\) is a band matrix with upper bandwidth
mu
and lower bandwidthml
, then the upper triangular factor \(U\) can have upper bandwidth as big assmu = MIN(N1,mu+ml)
. The lower triangular factor \(L\) has lower bandwidthml
.
The SUNLinSol_LapackBand module defines band implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_LapackBand
SUNLinSolInitialize_LapackBand
– this does nothing, since all consistency checks are performed at solver creation.SUNLinSolSetup_LapackBand
– this calls eitherDGBTRF
orSGBTRF
to perform the \(LU\) factorization.SUNLinSolSolve_LapackBand
– this calls eitherDGBTRS
orSGBTRS
to use the \(LU\) factors andpivots
array to perform the solve.SUNLinSolLastFlag_LapackBand
SUNLinSolSpace_LapackBand
– this only returns information for the storage within the solver object, i.e. storage forN
,last_flag
, andpivots
.SUNLinSolFree_LapackBand
10.12. The SUNLinSol_LapackDense Module
The SUNLinSol_LapackDense implementation of the SUNLinearSolver
class
is designed to be used with the corresponding SUNMATRIX_DENSE matrix type,
and one of the serial or sharedmemory N_Vector
implementations
(NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS).
10.12.1. SUNLinSol_LapackDense Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_lapackdense.h
. The installed module
library to link to is libsundials_sunlinsollapackdense
.lib
where .lib is typically .so
for shared libraries and
.a
for static libraries.
The module SUNLinSol_LapackDense provides the following additional usercallable constructor routine:

SUNLinearSolver SUNLinSol_LapackDense(N_Vector y, SUNMatrix A, SUNContext sunctx)
This function creates and allocates memory for a LAPACK dense
SUNLinearSolver
. Arguments:
y – vector used to determine the linear system size.
A – matrix used to assess compatibility.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
New SUNLinSol_LapackDense object, or
NULL
if eitherA
ory
are incompatible. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_DENSE matrix type and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.
10.12.2. SUNLinSol_LapackDense Description
The SUNLinSol_LapackDense module defines the
content field of a SUNLinearSolver
to be the following
structure:
struct _SUNLinearSolverContent_Dense {
sunindextype N;
sunindextype *pivots;
sunindextype last_flag;
};
These entries of the content field contain the following information:
N
 size of the linear system,pivots
 index array for partial pivoting in LU factorization,last_flag
 last error return flag from internal function evaluations.
The SUNLinSol_LapackDense module is a SUNLinearSolver
wrapper for
the LAPACK dense matrix factorization and solve routines, *GETRF
and *GETRS
, where *
is either D
or S
, depending on
whether SUNDIALS was configured to have sunrealtype
set to
double
or single
, respectively (see
§1.3 for details). In order to use the
SUNLinSol_LapackDense module it is assumed that LAPACK has been
installed on the system prior to installation of
SUNDIALS, and that SUNDIALS has been configured appropriately to
link with LAPACK (see
§1.2.4 for details).
We note that since there do not exist 128bit floatingpoint
factorization and solve routines in LAPACK, this interface cannot be
compiled when using extended
precision for sunrealtype
.
Similarly, since there do not exist 64bit integer LAPACK routines,
the SUNLinSol_LapackDense module also cannot be compiled when using
int64_t
for the sunindextype
.
This solver is constructed to perform the following operations:
The “setup” call performs an \(LU\) factorization with partial (row) pivoting (\(\mathcal O(N^3)\) cost), \(PA=LU\), where \(P\) is a permutation matrix, \(L\) is a lower triangular matrix with 1’s on the diagonal, and \(U\) is an upper triangular matrix. This factorization is stored inplace on the input SUNMATRIX_DENSE object \(A\), with pivoting information encoding \(P\) stored in the
pivots
array.The “solve” call performs pivoting and forward and backward substitution using the stored
pivots
array and the \(LU\) factors held in the SUNMATRIX_DENSE object (\(\mathcal O(N^2)\) cost).
The SUNLinSol_LapackDense module defines dense implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_LapackDense
SUNLinSolInitialize_LapackDense
– this does nothing, since all consistency checks are performed at solver creation.SUNLinSolSetup_LapackDense
– this calls eitherDGETRF
orSGETRF
to perform the \(LU\) factorization.SUNLinSolSolve_LapackDense
– this calls eitherDGETRS
orSGETRS
to use the \(LU\) factors andpivots
array to perform the solve.SUNLinSolLastFlag_LapackDense
SUNLinSolSpace_LapackDense
– this only returns information for the storage within the solver object, i.e. storage forN
,last_flag
, andpivots
.SUNLinSolFree_LapackDense
10.13. The SUNLinSol_MagmaDense Module
The SUNLinearSolver_MagmaDense implementation of the SUNLinearSolver
class is
designed to be used with the SUNMATRIX_MAGMADENSE matrix, and a GPUenabled
vector. The header file to include when using this module is
sunlinsol/sunlinsol_magmadense.h
. The installed library to link to is
libsundials_sunlinsolmagmadense.lib
where lib
is typically .so
for
shared libraries and .a
for static libraries.
Warning
The SUNLinearSolver_MagmaDense module is experimental and subject to change.
10.13.1. SUNLinearSolver_MagmaDense Description
The SUNLinearSolver_MagmaDense implementation provides an interface to the dense LU and dense batched LU methods in the MAGMA linear algebra library [126]. The batched LU methods are leveraged when solving block diagonal linear systems of the form
10.13.2. SUNLinearSolver_MagmaDense Functions
The SUNLinearSolver_MagmaDense module defines implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_MagmaDense
SUNLinSolInitialize_MagmaDense
SUNLinSolSetup_MagmaDense
SUNLinSolSolve_MagmaDense
SUNLinSolLastFlag_MagmaDense
SUNLinSolFree_MagmaDense
In addition, the module provides the following usercallable routines:

SUNLinearSolver SUNLinSol_MagmaDense(N_Vector y, SUNMatrix A, SUNContext sunctx)
This constructor function creates and allocates memory for a
SUNLinearSolver
object. Arguments:
y – a vector for checking compatibility with the solver.
A – a SUNMATRIX_MAGMADENSE matrix for checking compatibility with the solver.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either A or y are incompatible then this routine will returnNULL
. This routine analyzes the input matrix and vector to determine the linear system size and to assess compatibility with the solver.

SUNErrCode SUNLinSol_MagmaDense_SetAsync(SUNLinearSolver LS, sunbooleantype onoff)
This function can be used to toggle the linear solver between asynchronous and synchronous modes. In asynchronous mode (default), SUNLinearSolver operations are asynchronous with respect to the host. In synchronous mode, the host and GPU device are synchronized prior to the operation returning.
 Arguments:
LS – a SUNLinSol_MagmaDense object
onoff – 0 for synchronous mode or 1 for asynchronous mode (default 1)
 Return value:
10.13.3. SUNLinearSolver_MagmaDense Content
The SUNLinearSolver_MagmaDense module defines the object content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_MagmaDense {
int last_flag;
sunbooleantype async;
sunindextype N;
SUNMemory pivots;
SUNMemory pivotsarr;
SUNMemory dpivotsarr;
SUNMemory infoarr;
SUNMemory rhsarr;
SUNMemoryHelper memhelp;
magma_queue_t q;
};
10.14. The SUNLinSol_OneMklDense Module
The SUNLinearSolver_OneMklDense implementation of the SUNLinearSolver
class
interfaces to the direct linear solvers from the
Intel oneAPI Math Kernel Library (oneMKL)
for solving dense systems or blockdiagonal systems with dense blocks. This
linear solver is best paired with the SUNMatrix_OneMklDense matrix.
The header file to include when using this class is
sunlinsol/sunlinsol_onemkldense.h
. The installed library to link to is
libsundials_sunlinsolonemkldense.lib
where lib
is typically .so
for
shared libraries and .a
for static libraries.
Warning
The SUNLinearSolver_OneMklDense class is experimental and subject to change.
10.14.1. SUNLinearSolver_OneMklDense Functions
The SUNLinearSolver_OneMklDense class defines implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_OneMklDense
– returnsSUNLINEARSOLVER_ONEMKLDENSE
SUNLinSolInitialize_OneMklDense
SUNLinSolSetup_OneMklDense
SUNLinSolSolve_OneMklDense
SUNLinSolLastFlag_OneMklDense
SUNLinSolFree_OneMklDense
In addition, the class provides the following usercallable routines:

SUNLinearSolver SUNLinSol_OneMklDense(N_Vector y, SUNMatrix A, SUNContext sunctx)
This constructor function creates and allocates memory for a
SUNLinearSolver
object. Arguments:
y – a vector for checking compatibility with the solver.
A – a SUNMatrix_OneMklDense matrix for checking compatibility with the solver.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either A or y are incompatible then this routine will returnNULL
. This routine analyzes the input matrix and vector to determine the linear system size and to assess compatibility with the solver.
10.14.2. SUNLinearSolver_OneMklDense Usage Notes
Warning
The SUNLinearSolver_OneMklDense class only supports 64bit indexing, thus SUNDIALS must be built for 64bit indexing to use this class.
When using the SUNLinearSolver_OneMklDense class with a SUNDIALS package (e.g. CVODE), the queue given to the matrix is also used for the linear solver.
10.15. The SUNLinSol_PCG Module
The SUNLinSol_PCG implementation of the SUNLinearSolver
class performs
the PCG (Preconditioned Conjugate Gradient [65]) method; this
is an iterative linear solver that is designed to be compatible with
any N_Vector
implementation that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, and
N_VDestroy()
). Unlike the SPGMR and SPFGMR algorithms, PCG
requires a fixed amount of memory that does not increase with the
number of allowed iterations.
Unlike all of the other iterative linear solvers supplied with SUNDIALS, PCG should only be used on symmetric linear systems (e.g. mass matrix linear systems encountered in ARKODE). As a result, the explanation of the role of scaling and preconditioning matrices given in general must be modified in this scenario. The PCG algorithm solves a linear system \(Ax = b\) where \(A\) is a symmetric (\(A^T=A\)), realvalued matrix. Preconditioning is allowed, and is applied in a symmetric fashion on both the right and left. Scaling is also allowed and is applied symmetrically. We denote the preconditioner and scaling matrices as follows:
\(P\) is the preconditioner (assumed symmetric),
\(S\) is a diagonal matrix of scale factors.
The matrices \(A\) and \(P\) are not required explicitly; only routines
that provide \(A\) and \(P^{1}\) as operators are required. The diagonal
of the matrix \(S\) is held in a single N_Vector
, supplied by the user.
In this notation, PCG applies the underlying CG algorithm to the equivalent transformed system
where
The scaling matrix must be chosen so that the vectors \(SP^{1}b\) and \(S^{1}Px\) have dimensionless components.
The stopping test for the PCG iterations is on the L2 norm of the scaled preconditioned residual:
where \(\ v \_S = \sqrt{v^T S^T S v}\), with an input tolerance \(\delta\).
10.15.1. SUNLinSol_PCG Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_pcg.h
. The SUNLinSol_PCG module
is accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolpcg
module library.
The module SUNLinSol_PCG provides the following usercallable routines:

SUNLinearSolver SUNLinSol_PCG(N_Vector y, int pretype, int maxl, SUNContext sunctx)
This constructor function creates and allocates memory for a PCG
SUNLinearSolver
. Arguments:
y – a template vector.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
maxl – the maximum number of linear iterations to allow.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either y is incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with a consistent
N_Vector
implementation (i.e. that it supplies the requisite vector operations).A
maxl
argument that is \(\le0\) will result in the default value (5).Since the PCG algorithm is designed to only support symmetric preconditioning, then any of the
pretype
inputsSUN_PREC_LEFT
,SUN_PREC_RIGHT
, orSUN_PREC_BOTH
will result in use of the symmetric preconditioner; any other integer input will result in the default (no preconditioning). Although some SUNDIALS solvers are designed to only work with left preconditioning (IDA and IDAS) and others with only right preconditioning (KINSOL), PCG should only be used with these packages when the linear systems are known to be symmetric. Since the scaling of matrix rows and columns must be identical in a symmetric matrix, symmetric preconditioning should work appropriately even for packages designed with onesided preconditioning in mind.

SUNErrCode SUNLinSol_PCGSetPrecType(SUNLinearSolver S, int pretype)
This function updates the flag indicating use of preconditioning.
 Arguments:
S – SUNLinSol_PCG object to update.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
 Return value:
 Notes:
As above, any one of the input values,
SUN_PREC_LEFT
,SUN_PREC_RIGHT
, orSUN_PREC_BOTH
will enable preconditioning;SUN_PREC_NONE
disables preconditioning.

SUNErrCode SUNLinSol_PCGSetMaxl(SUNLinearSolver S, int maxl)
This function updates the number of linear solver iterations to allow.
 Arguments:
S – SUNLinSol_PCG object to update.
maxl – maximum number of linear iterations to allow. Any nonpositive input will result in the default value (5).
 Return value:
10.15.2. SUNLinSol_PCG Description
The SUNLinSol_PCG module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_PCG {
int maxl;
int pretype;
sunbooleantype zeroguess;
int numiters;
sunrealtype resnorm;
int last_flag;
SUNATimesFn ATimes;
void* ATData;
SUNPSetupFn Psetup;
SUNPSolveFn Psolve;
void* PData;
N_Vector s;
N_Vector r;
N_Vector p;
N_Vector z;
N_Vector Ap;
};
These entries of the content field contain the following information:
maxl
 number of PCG iterations to allow (default is 5),pretype
 flag for use of preconditioning (default is none),numiters
 number of iterations from the mostrecent solve,resnorm
 final linear residual norm from the mostrecent solve,last_flag
 last error return flag from an internal function,ATimes
 function pointer to perform \(Av\) product,ATData
 pointer to structure forATimes
,Psetup
 function pointer to preconditioner setup routine,Psolve
 function pointer to preconditioner solve routine,PData
 pointer to structure forPsetup
andPsolve
,s
 vector pointer for supplied scaling matrix (default isNULL
),r
 aN_Vector
which holds the preconditioned linear system residual,p, z, Ap
N_Vector
used for workspace by the PCG algorithm.
This solver is constructed to perform the following operations:
During construction all
N_Vector
solver data is allocated, with vectors cloned from a templateN_Vector
that is input, and default solver parameters are set.Userfacing “set” routines may be called to modify default solver parameters.
Additional “set” routines are called by the SUNDIALS solver that interfaces with SUNLinSol_PCG to supply the
ATimes
,PSetup
, andPsolve
function pointers ands
scaling vector.In the “initialize” call, the solver parameters are checked for validity.
In the “setup” call, any non
NULL
PSetup
function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the genericPSetup
function and the solverspecific routine (solversupplied or usersupplied).In the “solve” call the PCG iteration is performed. This will include scaling and preconditioning if those options have been supplied.
The SUNLinSol_PCG module defines implementations of all “iterative” linear solver operations listed in §10.1:
SUNLinSolGetType_PCG
SUNLinSolInitialize_PCG
SUNLinSolSetATimes_PCG
SUNLinSolSetPreconditioner_PCG
SUNLinSolSetScalingVectors_PCG
– since PCG only supports symmetric scaling, the secondN_Vector
argument to this function is ignored.SUNLinSolSetZeroGuess_PCG
– note the solver assumes a nonzero guess by default and the zero guess flag is reset toSUNFALSE
after each call toSUNLinSolSolve_PCG
.SUNLinSolSetup_PCG
SUNLinSolSolve_PCG
SUNLinSolNumIters_PCG
SUNLinSolResNorm_PCG
SUNLinSolResid_PCG
SUNLinSolLastFlag_PCG
SUNLinSolSpace_PCG
SUNLinSolFree_PCG
10.16. The SUNLinSol_SPBCGS Module
The SUNLinSol_SPBCGS implementation of the SUNLinearSolver
class performs
a Scaled, Preconditioned, BiConjugate Gradient, Stabilized [131] method;
this is an iterative linear solver that is designed to be compatible with any
N_Vector
implementation that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, N_VDiv()
, and
N_VDestroy()
). Unlike the SPGMR and SPFGMR algorithms,
SPBCGS requires a fixed amount of memory that does not increase with
the number of allowed iterations.
10.16.1. SUNLinSol_SPBCGS Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_spbcgs.h
. The SUNLinSol_SPBCGS module
is accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolspbcgs
module library.
The module SUNLinSol_SPBCGS provides the following usercallable routines:

SUNLinearSolver SUNLinSol_SPBCGS(N_Vector y, int pretype, int maxl, SUNContext sunctx)
This constructor function creates and allocates memory for a SPBCGS
SUNLinearSolver
. Arguments:
y – a template vector.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
maxl – the maximum number of linear iterations to allow.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either y is incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with a consistent
N_Vector
implementation (i.e. that it supplies the requisite vector operations).A
maxl
argument that is \(\le0\) will result in the default value (5).Some SUNDIALS solvers are designed to only work with left preconditioning (IDA and IDAS) and others with only right preconditioning (KINSOL). While it is possible to configure a SUNLinSol_SPBCGS object to use any of the preconditioning options with these solvers, this use mode is not supported and may result in inferior performance.
Note
With
SUN_PREC_RIGHT
orSUN_PREC_BOTH
the initial guess must be zero (useSUNLinSolSetZeroGuess()
to indicate the initial guess is zero).

SUNErrCode SUNLinSol_SPBCGSSetPrecType(SUNLinearSolver S, int pretype)
This function updates the flag indicating use of preconditioning.
 Arguments:
S – SUNLinSol_SPBCGS object to update.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
 Return value:

SUNErrCode SUNLinSol_SPBCGSSetMaxl(SUNLinearSolver S, int maxl)
This function updates the number of linear solver iterations to allow.
 Arguments:
S – SUNLinSol_SPBCGS object to update.
maxl – maximum number of linear iterations to allow. Any nonpositive input will result in the default value (5).
 Return value:
10.16.2. SUNLinSol_SPBCGS Description
The SUNLinSol_SPBCGS module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SPBCGS {
int maxl;
int pretype;
sunbooleantype zeroguess;
int numiters;
sunrealtype resnorm;
int last_flag;
SUNATimesFn ATimes;
void* ATData;
SUNPSetupFn Psetup;
SUNPSolveFn Psolve;
void* PData;
N_Vector s1;
N_Vector s2;
N_Vector r;
N_Vector r_star;
N_Vector p;
N_Vector q;
N_Vector u;
N_Vector Ap;
N_Vector vtemp;
};
These entries of the content field contain the following information:
maxl
 number of SPBCGS iterations to allow (default is 5),pretype
 flag for type of preconditioning to employ (default is none),numiters
 number of iterations from the mostrecent solve,resnorm
 final linear residual norm from the mostrecent solve,last_flag
 last error return flag from an internal function,ATimes
 function pointer to perform \(Av\) product,ATData
 pointer to structure forATimes
,Psetup
 function pointer to preconditioner setup routine,Psolve
 function pointer to preconditioner solve routine,PData
 pointer to structure forPsetup
andPsolve
,s1, s2
 vector pointers for supplied scaling matrices (default isNULL
),r
 aN_Vector
which holds the current scaled, preconditioned linear system residual,r_star
 aN_Vector
which holds the initial scaled, preconditioned linear system residual,p, q, u, Ap, vtemp
N_Vector
used for workspace by the SPBCGS algorithm.
This solver is constructed to perform the following operations:
During construction all
N_Vector
solver data is allocated, with vectors cloned from a templateN_Vector
that is input, and default solver parameters are set.Userfacing “set” routines may be called to modify default solver parameters.
Additional “set” routines are called by the SUNDIALS solver that interfaces with SUNLinSol_SPBCGS to supply the
ATimes
,PSetup
, andPsolve
function pointers ands1
ands2
scaling vectors.In the “initialize” call, the solver parameters are checked for validity.
In the “setup” call, any non
NULL
PSetup
function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the genericPSetup
function and the solverspecific routine (solversupplied or usersupplied).In the “solve” call the SPBCGS iteration is performed. This will include scaling and preconditioning if those options have been supplied.
The SUNLinSol_SPBCGS module defines implementations of all “iterative” linear solver operations listed in §10.1:
SUNLinSolGetType_SPBCGS
SUNLinSolInitialize_SPBCGS
SUNLinSolSetATimes_SPBCGS
SUNLinSolSetPreconditioner_SPBCGS
SUNLinSolSetScalingVectors_SPBCGS
SUNLinSolSetZeroGuess_SPBCGS
– note the solver assumes a nonzero guess by default and the zero guess flag is reset toSUNFALSE
after each call toSUNLinSolSolve_SPBCGS
.SUNLinSolSetup_SPBCGS
SUNLinSolSolve_SPBCGS
SUNLinSolNumIters_SPBCGS
SUNLinSolResNorm_SPBCGS
SUNLinSolResid_SPBCGS
SUNLinSolLastFlag_SPBCGS
SUNLinSolSpace_SPBCGS
SUNLinSolFree_SPBCGS
10.17. The SUNLinSol_SPFGMR Module
The SUNLinSol_SPFGMR implementation of the SUNLinearSolver
class performs
a Scaled, Preconditioned, Flexible, Generalized Minimum Residual [105]
method; this is an iterative linear solver that is designed to be compatible with
any N_Vector
implementation that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, N_VConst()
,
N_VDiv()
, and N_VDestroy()
). Unlike the other
Krylov iterative linear solvers supplied with SUNDIALS, FGMRES is
specifically designed to work with a changing preconditioner
(e.g. from an iterative method).
10.17.1. SUNLinSol_SPFGMR Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_spfgmr.h
. The SUNLinSol_SPFGMR module is
accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolspfgmr
module library.
The module SUNLinSol_SPFGMR provides the following usercallable routines:

SUNLinearSolver SUNLinSol_SPFGMR(N_Vector y, int pretype, int maxl, SUNContext sunctx)
This constructor function creates and allocates memory for a SPFGMR
SUNLinearSolver
. Arguments:
y – a template vector.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
maxl – the number of Krylov basis vectors to use.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either y is incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with a consistent
N_Vector
implementation (i.e. that it supplies the requisite vector operations).A
maxl
argument that is \(\le0\) will result in the default value (5).Since the FGMRES algorithm is designed to only support right preconditioning, then any of the
pretype
inputsSUN_PREC_LEFT
,SUN_PREC_RIGHT
, orSUN_PREC_BOTH
will result in use ofSUN_PREC_RIGHT
; any other integer input will result in the default (no preconditioning). We note that some SUNDIALS solvers are designed to only work with left preconditioning (IDA and IDAS). While it is possible to use a rightpreconditioned SUNLinSol_SPFGMR object for these packages, this use mode is not supported and may result in inferior performance.

SUNErrCode SUNLinSol_SPFGMRSetPrecType(SUNLinearSolver S, int pretype)
This function updates the flag indicating use of preconditioning.
 Arguments:
S – SUNLinSol_SPFGMR object to update.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
 Return value:
 Notes:
Since the FGMRES algorithm is designed to only support right preconditioning, then any of the
pretype
inputsSUN_PREC_LEFT
,SUN_PREC_RIGHT
, orSUN_PREC_BOTH
will result in use ofSUN_PREC_RIGHT
; any other integer input will result in the default (no preconditioning).

SUNErrCode SUNLinSol_SPFGMRSetGSType(SUNLinearSolver S, int gstype)
This function sets the type of GramSchmidt orthogonalization to use.
 Arguments:
S – SUNLinSol_SPFGMR object to update.
gstype – a flag indicating the type of orthogonalization to use:
SUN_MODIFIED_GS
SUN_CLASSICAL_GS
 Return value:

SUNErrCode SUNLinSol_SPFGMRSetMaxRestarts(SUNLinearSolver S, int maxrs)
This function sets the number of FGMRES restarts to allow.
 Arguments:
S – SUNLinSol_SPFGMR object to update.
maxrs – maximum number of restarts to allow. A negative input will result in the default of 0.
 Return value:
10.17.2. SUNLinSol_SPFGMR Description
The SUNLinSol_SPFGMR module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SPFGMR {
int maxl;
int pretype;
int gstype;
int max_restarts;
sunbooleantype zeroguess;
int numiters;
sunrealtype resnorm;
int last_flag;
SUNATimesFn ATimes;
void* ATData;
SUNPSetupFn Psetup;
SUNPSolveFn Psolve;
void* PData;
N_Vector s1;
N_Vector s2;
N_Vector *V;
N_Vector *Z;
sunrealtype **Hes;
sunrealtype *givens;
N_Vector xcor;
sunrealtype *yg;
N_Vector vtemp;
};
These entries of the content field contain the following information:
maxl
 number of FGMRES basis vectors to use (default is 5),pretype
 flag for use of preconditioning (default is none),gstype
 flag for type of GramSchmidt orthogonalization (default is modified GramSchmidt),max_restarts
 number of FGMRES restarts to allow (default is 0),numiters
 number of iterations from the mostrecent solve,resnorm
 final linear residual norm from the mostrecent solve,last_flag
 last error return flag from an internal function,ATimes
 function pointer to perform \(Av\) product,ATData
 pointer to structure forATimes
,Psetup
 function pointer to preconditioner setup routine,Psolve
 function pointer to preconditioner solve routine,PData
 pointer to structure forPsetup
andPsolve
,s1, s2
 vector pointers for supplied scaling matrices (default isNULL
),V
 the array of Krylov basis vectors \(v_1, \ldots, v_{\text{maxl}+1}\), stored inV[0], ..., V[maxl]
. Each \(v_i\) is a vector of typeN_Vector
,Z
 the array of preconditioned Krylov basis vectors \(z_1, \ldots, z_{\text{maxl}+1}\), stored inZ[0], ..., Z[maxl]
. Each \(z_i\) is a vector of typeN_Vector
,Hes
 the \((\text{maxl}+1)\times\text{maxl}\) Hessenberg matrix. It is stored rowwise so that the (i,j)th element is given byHes[i][j]
,givens
 a length \(2\,\text{maxl}\) array which represents the Givens rotation matrices that arise in the FGMRES algorithm. These matrices are \(F_0, F_1, \ldots, F_j\), where\[\begin{split}F_i = \begin{bmatrix} 1 & & & & & & & \\ & \ddots & & & & & & \\ & & 1 & & & & & \\ & & & c_i & s_i & & & \\ & & & s_i & c_i & & & \\ & & & & & 1 & & \\ & & & & & & \ddots & \\ & & & & & & & 1\end{bmatrix},\end{split}\]are represented in the
givens
vector asgivens[0]
\(= c_0\),givens[1]
\(= s_0\),givens[2]
\(= c_1\),givens[3]
\(= s_1\), \(\ldots\),givens[2j]
\(= c_j\),givens[2j+1]
\(= s_j\),xcor
 a vector which holds the scaled, preconditioned correction to the initial guess,yg
 a length \((\text{maxl}+1)\) array ofsunrealtype
values used to hold “short” vectors (e.g. \(y\) and \(g\)),vtemp
 temporary vector storage.
This solver is constructed to perform the following operations:
During construction, the
xcor
andvtemp
arrays are cloned from a templateN_Vector
that is input, and default solver parameters are set.Userfacing “set” routines may be called to modify default solver parameters.
Additional “set” routines are called by the SUNDIALS solver that interfaces with SUNLinSol_SPFGMR to supply the
ATimes
,PSetup
, andPsolve
function pointers ands1
ands2
scaling vectors.In the “initialize” call, the remaining solver data is allocated (
V
,Hes
,givens
, andyg
)In the “setup” call, any non
NULL
PSetup
function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the genericPSetup
function and the solverspecific routine (solversupplied or usersupplied).In the “solve” call, the FGMRES iteration is performed. This will include scaling, preconditioning, and restarts if those options have been supplied.
The SUNLinSol_SPFGMR module defines implementations of all “iterative” linear solver operations listed in §10.1:
SUNLinSolGetType_SPFGMR
SUNLinSolInitialize_SPFGMR
SUNLinSolSetATimes_SPFGMR
SUNLinSolSetPreconditioner_SPFGMR
SUNLinSolSetScalingVectors_SPFGMR
SUNLinSolSetZeroGuess_SPFGMR
– note the solver assumes a nonzero guess by default and the zero guess flag is reset toSUNFALSE
after each call toSUNLinSolSolve_SPFGMR
.SUNLinSolSetup_SPFGMR
SUNLinSolSolve_SPFGMR
SUNLinSolNumIters_SPFGMR
SUNLinSolResNorm_SPFGMR
SUNLinSolResid_SPFGMR
SUNLinSolLastFlag_SPFGMR
SUNLinSolSpace_SPFGMR
SUNLinSolFree_SPFGMR
10.18. The SUNLinSol_SPGMR Module
The SUNLinSol_SPGMR implementation of the SUNLinearSolver
class performs
a Scaled, Preconditioned, Generalized Minimum Residual [106] method;
this is an iterative linear solver that is designed to be compatible with any
N_Vector
implementation that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, N_VConst()
,
N_VDiv()
, and N_VDestroy()
).
10.18.1. SUNLinSol_SPGMR Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_spgmr.h
. The SUNinSol_SPGMR module
is accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolspgmr
module library.
The module SUNLinSol_SPGMR provides the following usercallable routines:

SUNLinearSolver SUNLinSol_SPGMR(N_Vector y, int pretype, int maxl, SUNContext sunctx)
This constructor function creates and allocates memory for a SPGMR
SUNLinearSolver
. Arguments:
y – a template vector.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
maxl – the number of Krylov basis vectors to use.
 Return value:
If successful, a
SUNLinearSolver
object. If either y is incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with a consistent
N_Vector
implementation (i.e. that it supplies the requisite vector operations).A
maxl
argument that is \(\le0\) will result in the default value (5).Some SUNDIALS solvers are designed to only work with left preconditioning (IDA and IDAS) and others with only right preconditioning (KINSOL). While it is possible to configure a SUNLinSol_SPGMR object to use any of the preconditioning options with these solvers, this use mode is not supported and may result in inferior performance.

SUNErrCode SUNLinSol_SPGMRSetPrecType(SUNLinearSolver S, int pretype)
This function updates the flag indicating use of preconditioning.
 Arguments:
S – SUNLinSol_SPGMR object to update.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
 Return value:

SUNErrCode SUNLinSol_SPGMRSetGSType(SUNLinearSolver S, int gstype)
This function sets the type of GramSchmidt orthogonalization to use.
 Arguments:
S – SUNLinSol_SPGMR object to update.
gstype – a flag indicating the type of orthogonalization to use:
SUN_MODIFIED_GS
SUN_CLASSICAL_GS
 Return value:

SUNErrCode SUNLinSol_SPGMRSetMaxRestarts(SUNLinearSolver S, int maxrs)
This function sets the number of GMRES restarts to allow.
 Arguments:
S – SUNLinSol_SPGMR object to update.
maxrs – maximum number of restarts to allow. A negative input will result in the default of 0.
 Return value:
10.18.2. SUNLinSol_SPGMR Description
The SUNLinSol_SPGMR module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SPGMR {
int maxl;
int pretype;
int gstype;
int max_restarts;
sunbooleantype zeroguess;
int numiters;
sunrealtype resnorm;
int last_flag;
SUNATimesFn ATimes;
void* ATData;
SUNPSetupFn Psetup;
SUNPSolveFn Psolve;
void* PData;
N_Vector s1;
N_Vector s2;
N_Vector *V;
sunrealtype **Hes;
sunrealtype *givens;
N_Vector xcor;
sunrealtype *yg;
N_Vector vtemp;
};
These entries of the content field contain the following information:
maxl
 number of GMRES basis vectors to use (default is 5),pretype
 flag for type of preconditioning to employ (default is none),gstype
 flag for type of GramSchmidt orthogonalization (default is modified GramSchmidt),max_restarts
 number of GMRES restarts to allow (default is 0),numiters
 number of iterations from the mostrecent solve,resnorm
 final linear residual norm from the mostrecent solve,last_flag
 last error return flag from an internal function,ATimes
 function pointer to perform \(Av\) product,ATData
 pointer to structure forATimes
,Psetup
 function pointer to preconditioner setup routine,Psolve
 function pointer to preconditioner solve routine,PData
 pointer to structure forPsetup
andPsolve
,s1, s2
 vector pointers for supplied scaling matrices (default isNULL
),V
 the array of Krylov basis vectors \(v_1, \ldots, v_{\text{maxl}+1}\), stored inV[0], ... V[maxl]
. Each \(v_i\) is a vector of typeN_Vector
,Hes
 the \((\text{maxl}+1)\times\text{maxl}\) Hessenberg matrix. It is stored rowwise so that the (i,j)th element is given byHes[i][j]
,givens
 a length \(2\,\text{maxl}\) array which represents the Givens rotation matrices that arise in the GMRES algorithm. These matrices are \(F_0, F_1, \ldots, F_j\), where\[\begin{split}F_i = \begin{bmatrix} 1 & & & & & & & \\ & \ddots & & & & & & \\ & & 1 & & & & & \\ & & & c_i & s_i & & & \\ & & & s_i & c_i & & & \\ & & & & & 1 & & \\ & & & & & & \ddots & \\ & & & & & & & 1\end{bmatrix},\end{split}\]are represented in the
givens
vector asgivens[0]
\(= c_0\),givens[1]
\(= s_0\),givens[2]
\(= c_1\),givens[3]
\(= s_1\), \(\ldots\),givens[2j]
\(= c_j\),givens[2j+1]
\(= s_j\),xcor
 a vector which holds the scaled, preconditioned correction to the initial guess,yg
 a length \((\text{maxl}+1)\) array ofsunrealtype
values used to hold “short” vectors (e.g. \(y\) and \(g\)),vtemp
 temporary vector storage.
This solver is constructed to perform the following operations:
During construction, the
xcor
andvtemp
arrays are cloned from a templateN_Vector
that is input, and default solver parameters are set.Userfacing “set” routines may be called to modify default solver parameters.
Additional “set” routines are called by the SUNDIALS solver that interfaces with SUNLinSol_SPGMR to supply the
ATimes
,PSetup
, andPsolve
function pointers ands1
ands2
scaling vectors.In the “initialize” call, the remaining solver data is allocated (
V
,Hes
,givens
, andyg
)In the “setup” call, any non
NULL
PSetup
function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the genericPSetup
function and the solverspecific routine (solversupplied or usersupplied).In the “solve” call, the GMRES iteration is performed. This will include scaling, preconditioning, and restarts if those options have been supplied.
The SUNLinSol_SPGMR module defines implementations of all “iterative” linear solver operations listed in §10.1:
SUNLinSolGetType_SPGMR
SUNLinSolInitialize_SPGMR
SUNLinSolSetATimes_SPGMR
SUNLinSolSetPreconditioner_SPGMR
SUNLinSolSetScalingVectors_SPGMR
SUNLinSolSetZeroGuess_SPGMR
– note the solver assumes a nonzero guess by default and the zero guess flag is reset toSUNFALSE
after each call toSUNLinSolSolve_SPGMR
.SUNLinSolSetup_SPGMR
SUNLinSolSolve_SPGMR
SUNLinSolNumIters_SPGMR
SUNLinSolResNorm_SPGMR
SUNLinSolResid_SPGMR
SUNLinSolLastFlag_SPGMR
SUNLinSolSpace_SPGMR
SUNLinSolFree_SPGMR
10.19. The SUNLinSol_SPTFQMR Module
The SUNLinSol_SPTFQMR implementation of the SUNLinearSolver
class performs
a Scaled, Preconditioned, TransposeFree QuasiMinimum Residual [56]
method; this is an iterative linear solver that is designed to be compatible
with any N_Vector
implementation that supports a minimal subset of operations
(N_VClone()
, N_VDotProd()
, N_VScale()
,
N_VLinearSum()
, N_VProd()
, N_VConst()
,
N_VDiv()
, and N_VDestroy()
). Unlike the SPGMR and
SPFGMR algorithms, SPTFQMR requires a fixed amount of memory that does
not increase with the number of allowed iterations.
10.19.1. SUNLinSol_SPTFQMR Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_sptfqmr.h
. The SUNLinSol_SPTFQMR module
is accessible from all SUNDIALS solvers without
linking to the libsundials_sunlinsolsptfqmr
module library.
The module SUNLinSol_SPTFQMR provides the following usercallable routines:

SUNLinearSolver SUNLinSol_SPTFQMR(N_Vector y, int pretype, int maxl, SUNContext sunctx)
This constructor function creates and allocates memory for a SPTFQMR
SUNLinearSolver
. Arguments:
y – a template vector.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
maxl – the number of Krylov basis vectors to use.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either y is incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with a consistent
N_Vector
implementation (i.e. that it supplies the requisite vector operations).A
maxl
argument that is \(\le0\) will result in the default value (5).Some SUNDIALS solvers are designed to only work with left preconditioning (IDA and IDAS) and others with only right preconditioning (KINSOL). While it is possible to configure a SUNLinSol_SPTFQMR object to use any of the preconditioning options with these solvers, this use mode is not supported and may result in inferior performance.
Note
With
SUN_PREC_RIGHT
orSUN_PREC_BOTH
the initial guess must be zero (useSUNLinSolSetZeroGuess()
to indicate the initial guess is zero).

SUNErrCode SUNLinSol_SPTFQMRSetPrecType(SUNLinearSolver S, int pretype)
This function updates the flag indicating use of preconditioning.
 Arguments:
S – SUNLinSol_SPGMR object to update.
pretype – a flag indicating the type of preconditioning to use:
SUN_PREC_NONE
SUN_PREC_LEFT
SUN_PREC_RIGHT
SUN_PREC_BOTH
 Return value:

SUNErrCode SUNLinSol_SPTFQMRSetMaxl(SUNLinearSolver S, int maxl)
This function updates the number of linear solver iterations to allow.
 Arguments:
S – SUNLinSol_SPTFQMR object to update.
maxl – maximum number of linear iterations to allow. Any nonpositive input will result in the default value (5).
 Return value:
10.19.2. SUNLinSol_SPTFQMR Description
The SUNLinSol_SPTFQMR module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SPTFQMR {
int maxl;
int pretype;
sunbooleantype zeroguess;
int numiters;
sunrealtype resnorm;
int last_flag;
SUNATimesFn ATimes;
void* ATData;
SUNPSetupFn Psetup;
SUNPSolveFn Psolve;
void* PData;
N_Vector s1;
N_Vector s2;
N_Vector r_star;
N_Vector q;
N_Vector d;
N_Vector v;
N_Vector p;
N_Vector *r;
N_Vector u;
N_Vector vtemp1;
N_Vector vtemp2;
N_Vector vtemp3;
};
These entries of the content field contain the following information:
maxl
 number of TFQMR iterations to allow (default is 5),pretype
 flag for type of preconditioning to employ (default is none),numiters
 number of iterations from the mostrecent solve,resnorm
 final linear residual norm from the mostrecent solve,last_flag
 last error return flag from an internal function,ATimes
 function pointer to perform \(Av\) product,ATData
 pointer to structure forATimes
,Psetup
 function pointer to preconditioner setup routine,Psolve
 function pointer to preconditioner solve routine,PData
 pointer to structure forPsetup
andPsolve
,s1, s2
 vector pointers for supplied scaling matrices (default isNULL
),r_star
 aN_Vector
which holds the initial scaled, preconditioned linear system residual,q, d, v, p, u
N_Vector
used for workspace by the SPTFQMR algorithm,r
 array of twoN_Vector
used for workspace within the SPTFQMR algorithm,vtemp1, vtemp2, vtemp3
 temporary vector storage.
This solver is constructed to perform the following operations:
During construction all
N_Vector
solver data is allocated, with vectors cloned from a templateN_Vector
that is input, and default solver parameters are set.Userfacing “set” routines may be called to modify default solver parameters.
Additional “set” routines are called by the SUNDIALS solver that interfaces with SUNLinSol_SPTFQMR to supply the
ATimes
,PSetup
, andPsolve
function pointers ands1
ands2
scaling vectors.In the “initialize” call, the solver parameters are checked for validity.
In the “setup” call, any non
NULL
PSetup
function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the genericPSetup
function and the solverspecific routine (solversupplied or usersupplied).In the “solve” call the TFQMR iteration is performed. This will include scaling and preconditioning if those options have been supplied.
The SUNLinSol_SPTFQMR module defines implementations of all “iterative” linear solver operations listed in §10.1:
SUNLinSolGetType_SPTFQMR
SUNLinSolInitialize_SPTFQMR
SUNLinSolSetATimes_SPTFQMR
SUNLinSolSetPreconditioner_SPTFQMR
SUNLinSolSetScalingVectors_SPTFQMR
SUNLinSolSetZeroGuess_SPTFQMR
– note the solver assumes a nonzero guess by default and the zero guess flag is reset toSUNFALSE
after each call toSUNLinSolSolve_SPTFQMR
.SUNLinSolSetup_SPTFQMR
SUNLinSolSolve_SPTFQMR
SUNLinSolNumIters_SPTFQMR
SUNLinSolResNorm_SPTFQMR
SUNLinSolResid_SPTFQMR
SUNLinSolLastFlag_SPTFQMR
SUNLinSolSpace_SPTFQMR
SUNLinSolFree_SPTFQMR
10.20. The SUNLinSol_SuperLUDIST Module
The SUNLinsol_SuperLUDIST implementation of the SUNLinearSolver
class interfaces
with the SuperLU_DIST library. This is designed to be used with the
SUNMatrix_SLUNRloc SUNMatrix
, and one of the serial, threaded or parallel
N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP, NVECTOR_PTHREADS,
NVECTOR_PARALLEL, NVECTOR_PARHYP).
10.20.1. SUNLinSol_SuperLUDIST Usage
The header file to be included when using this module
is sunlinsol/sunlinsol_superludist.h
. The installed module
library to link to is libsundials_sunlinsolsuperludist
.lib
where .lib is typically .so
for shared libraries and
.a
for static libraries.
The module SUNLinSol_SuperLUDIST provides the following usercallable routines:
Warning
Starting with SuperLU_DIST version 6.3.0, some structures were
renamed to have a prefix for the floating point type. The double precision API
functions have the prefix ‘d’. To maintain backwards compatibility with the
unprefixed types, SUNDIALS provides macros to these SuperLU_DIST types with an
‘x’ prefix that expand to the correct prefix. E.g., the SUNDIALS macro
xLUstruct_t
expands to dLUstruct_t
or LUstruct_t
based on the
SuperLU_DIST version.

SUNLinearSolver SUNLinSol_SuperLUDIST(N_Vector y, SuperMatrix *A, gridinfo_t *grid, xLUstruct_t *lu, xScalePermstruct_t *scaleperm, xSOLVEstruct_t *solve, SuperLUStat_t *stat, superlu_dist_options_t *options, SUNContext sunctx)
This constructor function creates and allocates memory for a SUNLinSol_SuperLUDIST object.
 Arguments:
y – a template vector.
A – a template matrix
grid, lu, scaleperm, solve, stat, options – SuperLU_DIST object pointers.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object; otherwise this routine will returnNULL
. Notes:
This routine analyzes the input matrix and vector to determine the linear system size and to assess the compatibility with the SuperLU_DIST library.
This routine will perform consistency checks to ensure that it is called with consistent N_Vector and SUNMatrix implementations. These are currently limited to the SUNMatrix_SLUNRloc matrix type and the NVECTOR_SERIAL, NVECTOR_OPENMP, NVECTOR_PTHREADS, NVECTOR_PARALLEL, and NVECTOR_PARHYP vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.
The
grid
,lu
,scaleperm
,solve
, andoptions
arguments are not checked and are passed directly to SuperLU_DIST routines.Some struct members of the
options
argument are modified internally by the SUNLinSol_SuperLUDIST solver. Specifically, the memberFact
is modified in the setup and solve routines.

sunrealtype SUNLinSol_SuperLUDIST_GetBerr(SUNLinearSolver LS)
This function returns the componentwise relative backward error of the computed solution. It takes one argument, the
SUNLinearSolver
object. The return type issunrealtype
.

gridinfo_t *SUNLinSol_SuperLUDIST_GetGridinfo(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that contains the 2D process grid. It takes one argument, the
SUNLinearSolver
object.

xLUstruct_t *SUNLinSol_SuperLUDIST_GetLUstruct(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that contains the distributed
L
andU
structures. It takes one argument, theSUNLinearSolver
object.

superlu_dist_options_t *SUNLinSol_SuperLUDIST_GetSuperLUOptions(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that contains the options which control how the linear system is factorized and solved. It takes one argument, the
SUNLinearSolver
object.

xScalePermstruct_t *SUNLinSol_SuperLUDIST_GetScalePermstruct(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that contains the vectors that describe the transformations done to the matrix
A
. It takes one argument, theSUNLinearSolver
object.

xSOLVEstruct_t *SUNLinSol_SuperLUDIST_GetSOLVEstruct(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that contains information for communication during the solution phase. It takes one argument the
SUNLinearSolver
object.

SuperLUStat_t *SUNLinSol_SuperLUDIST_GetSuperLUStat(SUNLinearSolver LS)
This function returns a pointer to the SuperLU_DIST structure that stores information about runtime and flop count. It takes one argument, the
SUNLinearSolver
object.
10.20.2. SUNLinSol_SuperLUDIST Description
The SUNLinSol_SuperLUDIST module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SuperLUDIST {
sunbooleantype first_factorize;
int last_flag;
sunrealtype berr;
gridinfo_t *grid;
xLUstruct_t *lu;
superlu_dist_options_t *options;
xScalePermstruct_t *scaleperm;
xSOLVEstruct_t *solve;
SuperLUStat_t *stat;
sunindextype N;
};
These entries of the content field contain the following information:
first_factorize
– flag indicating whether the factorization has ever been performed,last_flag
– last error return flag from internal function evaluations,berr
– the componentwise relative backward error of the computed solution,grid
– pointer to the SuperLU_DIST structure that strores the 2D process gridlu
– pointer to the SuperLU_DIST structure that stores the distributedL
andU
factors,scaleperm
– pointer to the SuperLU_DIST structure that stores vectors describing the transformations done to the matrixA
,options
– pointer to the SuperLU_DIST stucture which contains options that control how the linear system is factorized and solved,solve
– pointer to the SuperLU_DIST solve structure,stat
– pointer to the SuperLU_DIST structure that stores information about runtime and flop count,N
– the number of equations in the system.
The SUNLinSol_SuperLUDIST module is a SUNLinearSolver adapter for the SuperLU_DIST sparse matrix factorization and solver library written by X. Sherry Li and collaborators [8, 58, 90, 91]. The package uses a SPMD parallel programming model and multithreading to enhance efficiency in distributedmemory parallel environments with multicore nodes and possibly GPU accelerators. It uses MPI for communication, OpenMP for threading, and CUDA for GPU support. In order to use the SUNLinSol_SuperLUDIST interface to SuperLU_DIST, it is assumed that SuperLU_DIST has been installed on the system prior to installation of SUNDIALS, and that SUNDIALS has been configured appropriately to link with SuperLU_DIST (see §1.2.4 for details). Additionally, the wrapper only supports doubleprecision calculations, and therefore cannot be compiled if SUNDIALS is configured to use single or extended precision. Moreover, since the SuperLU_DIST library may be installed to support either 32bit or 64bit integers, it is assumed that the SuperLU_DIST library is installed using the same integer size as SUNDIALS.
The SuperLU_DIST library provides many options to control how a linear
system will be factorized and solved. These options may be set by a user
on an instance of the superlu_dist_options_t
struct, and then it may be provided
as an argument to the SUNLinSol_SuperLUDIST constructor. The SUNLinSol_SuperLUDIST
module will respect all options set except for Fact
– this option is
necessarily modified by the SUNLinSol_SuperLUDIST module in the setup and solve routines.
Since the linear systems that arise within the context of SUNDIALS calculations will typically have identical sparsity patterns, the SUNLinSol_SuperLUDIST module is constructed to perform the following operations:
The first time that the “setup” routine is called, it sets the SuperLU_DIST option
Fact
toDOFACT
so that a subsequent call to the “solve” routine will perform a symbolic factorization, followed by an initial numerical factorization before continuing to solve the system.On subsequent calls to the “setup” routine, it sets the SuperLU_DIST option
Fact
toSamePattern
so that a subsequent call to “solve” will perform factorization assuming the same sparsity pattern as prior, i.e. it will reuse the column permutation vector.If “setup” is called prior to the “solve” routine, then the “solve” routine will perform a symbolic factorization, followed by an initial numerical factorization before continuing to the sparse triangular solves, and, potentially, iterative refinement. If “setup” is not called prior, “solve” will skip to the triangular solve step. We note that in this solve SuperLU_DIST operates on the native data arrays for the righthand side and solution vectors, without requiring costly data copies.
The SUNLinSol_SuperLUDIST module defines implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_SuperLUDIST
SUNLinSolInitialize_SuperLUDIST
– this sets thefirst_factorize
flag to 1 and resets the internal SuperLU_DIST statistics variables.SUNLinSolSetup_SuperLUDIST
– this sets the appropriate SuperLU_DIST options so that a subsequent solve will perform a symbolic and numerical factorization before proceeding with the triangular solvesSUNLinSolSolve_SuperLUDIST
– this calls the SuperLU_DIST solve routine to perform factorization (if the setup routine was called prior) and then use the $LU$ factors to solve the linear system.SUNLinSolLastFlag_SuperLUDIST
SUNLinSolSpace_SuperLUDIST
– this only returns information for the storage within the solver interface, i.e. storage for the integerslast_flag
andfirst_factorize
. For additional space requirements, see the SuperLU_DIST documentation.SUNLinSolFree_SuperLUDIST
10.21. The SUNLinSol_SuperLUMT Module
The SUNLinSol_SuperLUMT implementation of the SUNLinearSolver
class
interfaces with the SuperLU_MT library. This is designed to be used
with the corresponding SUNMATRIX_SPARSE matrix type, and one of the
serial or sharedmemory N_Vector
implementations (NVECTOR_SERIAL,
NVECTOR_OPENMP, or NVECTOR_PTHREADS). While these are compatible, it
is not recommended to use a threaded vector module with
SUNLinSol_SuperLUMT unless it is the NVECTOR_OPENMP module and the
SuperLU_MT library has also been compiled with OpenMP.
10.21.1. SUNLinSol_SuperLUMT Usage
The header file to be included when using this module
is sunlinsol/sunlinsol.SuperLUMT.h
. The installed module
library to link to is libsundials_sunlinsolsuperlumt
.lib
where .lib is typically .so
for shared libraries and
.a
for static libraries.
The module SUNLinSol_SuperLUMT provides the following usercallable routines:

SUNLinearSolver SUNLinSol_SuperLUMT(N_Vector y, SUNMatrix A, int num_threads, SUNContext sunctx)
This constructor function creates and allocates memory for a SUNLinSol_SuperLUMT object.
 Arguments:
y – a template vector.
A – a template matrix
num_threads – desired number of threads (OpenMP or Pthreads, depending on how SuperLU_MT was installed) to use during the factorization steps.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object; otherwise this routine will returnNULL
. Notes:
This routine analyzes the input matrix and vector to determine the linear system size and to assess compatibility with the SuperLU_MT library.
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_SPARSE matrix type (using either CSR or CSC storage formats) and the NVECTOR_SERIAL, NVECTOR_OPENMP, and NVECTOR_PTHREADS vector types. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.The
num_threads
argument is not checked and is passed directly to SuperLU_MT routines.

SUNErrCode SUNLinSol_SuperLUMTSetOrdering(SUNLinearSolver S, int ordering_choice)
This function sets the ordering used by SuperLU_MT for reducing fill in the linear solve.
 Arguments:
S – the SUNLinSol_SuperLUMT object to update.
ordering_choice:
natural ordering
minimal degree ordering on \(A^TA\)
minimal degree ordering on \(A^T+A\)
COLAMD ordering for unsymmetric matrices
The default is 3 for COLAMD.
 Return value:
10.21.2. SUNLinSol_SuperLUMT Description
The SUNLinSol_SuperLUMT module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_SuperLUMT {
int last_flag;
int first_factorize;
SuperMatrix *A, *AC, *L, *U, *B;
Gstat_t *Gstat;
sunindextype *perm_r, *perm_c;
sunindextype N;
int num_threads;
sunrealtype diag_pivot_thresh;
int ordering;
superlumt_options_t *options;
};
These entries of the content field contain the following information:
last_flag
 last error return flag from internal function evaluations,first_factorize
 flag indicating whether the factorization has ever been performed,A, AC, L, U, B
SuperMatrix
pointers used in solve,Gstat
GStat_t
object used in solve,perm_r, perm_c
 permutation arrays used in solve,N
 size of the linear system,num_threads
 number of OpenMP/Pthreads threads to use,diag_pivot_thresh
 threshold on diagonal pivoting,ordering
 flag for which reordering algorithm to use,options
 pointer to SuperLU_MT options structure.
The SUNLinSol_SuperLUMT module is a SUNLinearSolver
wrapper for
the SuperLU_MT sparse matrix factorization and solver library
written by X. Sherry Li and collaborators
[9, 42, 89]. The
package performs matrix factorization using threads to enhance
efficiency in shared memory parallel environments. It should be noted
that threads are only used in the factorization step. In
order to use the SUNLinSol_SuperLUMT interface to SuperLU_MT, it is
assumed that SuperLU_MT has been installed on the system prior to
installation of SUNDIALS, and that SUNDIALS has been configured
appropriately to link with SuperLU_MT (see
§1.2.4 for details).
Additionally, this wrapper only supports single and
doubleprecision calculations, and therefore cannot be compiled if
SUNDIALS is configured to have sunrealtype
set to extended
(see §1.3 for details). Moreover,
since the SuperLU_MT library may be installed to support either 32bit
or 64bit integers, it is assumed that the SuperLU_MT library is
installed using the same integer precision as the SUNDIALS
sunindextype
option.
The SuperLU_MT library has a symbolic factorization routine that computes the permutation of the linear system matrix to reduce fillin on subsequent \(LU\) factorizations (using COLAMD, minimal degree ordering on \(A^T*A\), minimal degree ordering on \(A^T+A\), or natural ordering). Of these ordering choices, the default value in the SUNLinSol_SuperLUMT module is the COLAMD ordering.
Since the linear systems that arise within the context of SUNDIALS calculations will typically have identical sparsity patterns, the SUNLinSol_SuperLUMT module is constructed to perform the following operations:
The first time that the “setup” routine is called, it performs the symbolic factorization, followed by an initial numerical factorization.
On subsequent calls to the “setup” routine, it skips the symbolic factorization, and only refactors the input matrix.
The “solve” call performs pivoting and forward and backward substitution using the stored SuperLU_MT data structures. We note that in this solve SuperLU_MT operates on the native data arrays for the righthand side and solution vectors, without requiring costly data copies.
The SUNLinSol_SuperLUMT module defines implementations of all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_SuperLUMT
SUNLinSolInitialize_SuperLUMT
– this sets thefirst_factorize
flag to 1 and resets the internal SuperLU_MT statistics variables.SUNLinSolSetup_SuperLUMT
– this performs either a \(LU\) factorization or refactorization of the input matrix.SUNLinSolSolve_SuperLUMT
– this calls the appropriate SuperLU_MT solve routine to utilize the \(LU\) factors to solve the linear system.SUNLinSolLastFlag_SuperLUMT
SUNLinSolSpace_SuperLUMT
– this only returns information for the storage within the solver interface, i.e. storage for the integerslast_flag
andfirst_factorize
. For additional space requirements, see the SuperLU_MT documentation.SUNLinSolFree_SuperLUMT
10.22. The SUNLinSol_cuSolverSp_batchQR Module
The SUNLinSol_cuSolverSp_batchQR implementation of the SUNLinearSolver
class
is designed to be used with the SUNMATRIX_CUSPARSE matrix, and the NVECTOR_CUDA vector.
The header file to include when using this module is
sunlinsol/sunlinsol_cusolversp_batchqr.h
. The installed library to link to
is libsundials_sunlinsolcusolversp.lib
where .lib
is typically
.so
for shared libraries and .a
for static libraries.
Warning
The SUNLinearSolver_cuSolverSp_batchQR module is experimental and subject to change.
10.22.1. SUNLinSol_cuSolverSp_batchQR description
The SUNLinearSolver_cuSolverSp_batchQR implementation provides an interface to the batched sparse QR factorization method provided by the NVIDIA cuSOLVER library [6]. The module is designed for solving block diagonal linear systems of the form
where all block matrices \(\mathbf{A_j}\) share the same sparsity pattern. The matrix
must be the SUNMatrix.cuSparse
.
10.22.2. SUNLinSol_cuSolverSp_batchQR functions
The SUNLinearSolver_cuSolverSp_batchQR
module defines implementations of
all “direct” linear solver operations listed in §10.1:
SUNLinSolGetType_cuSolverSp_batchQR
SUNLinSolInitialize_cuSolverSp_batchQR
– this sets thefirst_factorize
flag to 1SUNLinSolSetup_cuSolverSp_batchQR
– this always copies the relevant SUNMATRIX_SPARSE data to the GPU; if this is the first setup it will perform symbolic analysis on the systemSUNLinSolSolve_cuSolverSp_batchQR
– this calls thecusolverSpXcsrqrsvBatched
routine to perform factorizationSUNLinSolLastFlag_cuSolverSp_batchQR
SUNLinSolFree_cuSolverSp_batchQR
In addition, the module provides the following usercallable routines:

SUNLinearSolver SUNLinSol_cuSolverSp_batchQR(N_Vector y, SUNMatrix A, cusolverHandle_t cusol, SUNContext sunctx)
The function
SUNLinSol_cuSolverSp_batchQR
creates and allocates memory for a SUNLinearSolver object. Arguments:
y – a vector for checking compatibility with the solver.
A – a SUNMATRIX_cuSparse matrix for checking compatibility with the solver.
cusol – cuSolverSp object to use.
sunctx – the
SUNContext
object (see §1.4)
 Return value:
If successful, a
SUNLinearSolver
object. If either A or y are incompatible then this routine will returnNULL
. Notes:
This routine will perform consistency checks to ensure that it is called with consistent
N_Vector
andSUNMatrix
implementations. These are currently limited to the SUNMATRIX_CUSPARSE matrix type and the NVECTOR_CUDA vector type. Since the SUNMATRIX_CUSPARSE matrix type is only compatible with the NVECTOR_CUDA the restriction is also in place for the linear solver. As additional compatible matrix and vector implementations are added to SUNDIALS, these will be included within this compatibility check.

void SUNLinSol_cuSolverSp_batchQR_GetDescription(SUNLinearSolver LS, char **desc)
The function
SUNLinSol_cuSolverSp_batchQR_GetDescription
accesses the string description of the object (empty by default).

void SUNLinSol_cuSolverSp_batchQR_SetDescription(SUNLinearSolver LS, const char *desc)
The function
SUNLinSol_cuSolverSp_batchQR_SetDescription
sets the string description of the object (empty by default).

void SUNLinSol_cuSolverSp_batchQR_GetDeviceSpace(SUNLinearSolver S, size_t *cuSolverInternal, size_t *cuSolverWorkspace)
The function
SUNLinSol_cuSolverSp_batchQR_GetDeviceSpace
returns the cuSOLVER batch QR method internal buffer size, in bytes, in the argumentcuSolverInternal
and the cuSOLVER batch QR workspace buffer size, in bytes, in the agrumentcuSolverWorkspace
. The size of the internal buffer is proportional to the number of matrix blocks while the size of the workspace is almost independent of the number of blocks.
10.22.3. SUNLinSol_cuSolverSp_batchQR content
The SUNLinSol_cuSolverSp_batchQR module defines the content field of a
SUNLinearSolver
to be the following structure:
struct _SUNLinearSolverContent_cuSolverSp_batchQR {
int last_flag; /* last return flag */
sunbooleantype first_factorize; /* is this the first factorization? */
size_t internal_size; /* size of cusolver buffer for Q and R */
size_t workspace_size; /* size of cusolver memory for factorization */
cusolverSpHandle_t cusolver_handle; /* cuSolverSp context */
csrqrInfo_t info; /* opaque cusolver data structure */
void* workspace; /* memory block used by cusolver */
const char* desc; /* description of this linear solver */
};
10.23. The SUNLINEARSOLVER_GINKGO Module
Added in version 6.4.0.
The SUNLINEARSOLVER_GINKGO implementation of the SUNLinearSolver
API provides an
interface to the linear solvers from the Ginkgo linear algebra library [12].
Since Ginkgo is a modern C++ library, SUNLINEARSOLVER_GINKGO is also written in
modern C++ (specifically, C++14). Unlike most other SUNDIALS modules, it is
a header only library. To use the SUNLINEARSOLVER_GINKGO SUNLinearSolver
, users will
need to include sunlinsol/sunlinsol_ginkgo.hpp
. The module is meant to be used with
the SUNMATRIX_GINKGO module described in §9.16.
Instructions on building SUNDIALS with Ginkgo enabled are given
in §1.2.4. For instructions on
building and using Ginkgo itself, refer to the
Ginkgo website and documentation.
Note
It is assumed that users of this module are aware of how to use Ginkgo. This module does not try to encapsulate Ginkgo linear solvers, rather it provides a lightweight iteroperability layer between Ginkgo and SUNDIALS. Most, if not all, of the Ginkgo linear solver should work with this interface.
10.23.1. Using SUNLINEARSOLVER_GINKGO
After choosing a compatible N_Vector
(see §9.16.1) and creating a Ginkgoenabled SUNMatrix
(see
§9.16) to use the SUNLINEARSOLVER_GINKGO module, we first create a Ginkgo
stopping criteria object. Importantly, the sundials::ginkgo::DefaultStop
class provided
by SUNDIALS implements a stopping critierion that matches the default SUNDIALS stopping critierion.
Namely, it checks if the max iterations (5 by default) were reached or if the absolute residual
norm was below a speicified tolerance. The critierion can be created just like any other
Ginkgo stopping criteria:
auto crit{sundials::ginkgo::DefaultStop::build().with_max_iters(max_iters).on(gko_exec)};
Warning
It is highly recommended to employ this criterion when using Ginkgo solvers with SUNDIALS, but it is optional. However, to use the Ginkgo multigrid or cbgmres linear solvers, different Ginkgo criterion must be used.
Once we have created our stopping critierion, we create a Ginkgo solver factory object and
wrap it in a sundials::ginkgo::LinearSolver
object. In this example, we create
a Ginkgo conjugate gradient solver:
using GkoMatrixType = gko::matrix::Csr<sunrealtype, sunindextype>;
using GkoSolverType = gko::solver::Cg<sunrealtype>;
auto gko_solver_factory = gko::share(
GkoSolverType::build().with_criteria(std::move(crit)).on(gko_exec));
sundials::ginkgo::LinearSolver<GkoSolverType, GkoMatrixType> LS{
gko_solver_factory, sunctx};
Finally, we can pass the instance of sundials::ginkgo::LinearSolver
to any function
expecting a SUNLinearSolver
object through the implicit conversion operator or explicit conversion function.
// Attach linear solver and matrix to CVODE.
//
// Implicit conversion from sundials::ginkgo::LinearSolver<GkoSolverType, GkoMatrixType>
// to a SUNLinearSolver object is done.
//
// For details about creating A see the SUNMATRIX_GINKGO module.
CVodeSetLinearSolver(cvode_mem, LS, A);
// Alternatively with explicit conversion of LS to a SUNLinearSolver
// and A to a SUNMatrix:
CVodeSetLinearSolver(cvode_mem, LS>Convert(), A>Convert());
Warning
SUNLinSolFree()
should never be called on a SUNLinearSolver
that was created via conversion
from a sundials::ginkgo::LinearSolver
. Doing so may result in a double free.
10.23.2. SUNLINEARSOLVER_GINKGO API
In this section we list the public API of the sundials::ginkgo::LinearSolver
class.

template<class GkoSolverType, class GkoMatrixType>
class sundials::ginkgo::LinearSolver : public sundials::ConvertibleTo<SUNLinearSolver> 
LinearSolver() = default;
Default constructor  means the solver must be moved to.
Constructs a new LinearSolver from a Ginkgo solver factory.
 Parameters:
gko_solver_factory – The Ginkgo solver factory (typically gko::matrix::<type>::Factory`)
sunctx – The SUNDIALS simulation context (
SUNContext
)

LinearSolver(LinearSolver &&that_solver) noexcept
Move constructor.

LinearSolver &operator=(LinearSolver &&rhs)
Move assignment.

~LinearSolver() override = default
Default destructor.

operator SUNLinearSolver() override
Implicit conversion to a
SUNLinearSolver
.

operator SUNLinearSolver() const override
Implicit conversion to a
SUNLinearSolver
.

SUNLinearSolver Convert() override
Explicit conversion to a
SUNLinearSolver
.

SUNLinearSolver Convert() const override
Explicit conversion to a
SUNLinearSolver
.

std::shared_ptr<const gko::Executor> GkoExec() const
Get the
gko::Executor
associated with the Ginkgo solver.

std::shared_ptr<typename GkoSolverType::Factory> GkoFactory()
Get the underlying Ginkgo solver factory.

GkoSolverType *GkoSolver()
Get the underlying Ginkgo solver.
Note
This will be nullptr until the linear solver setup phase.

int NumIters() const
Get the number of linear solver iterations in the most recent solve.

sunrealtype ResNorm() const
Get the residual norm of the solution at the end of the last solve.
The type of residual norm depends on the Ginkgo stopping criteria used with the solver. With the
DefaultStop
criteria this would be the absolute residual 2norm.

GkoSolverType *Setup(Matrix<GkoMatrixType> *A)
Setup the linear system.
 Parameters:
A – the linear system matrix
 Returns:
Pointer to the Ginkgo solver generated from the factory

gko::LinOp *Solve(N_Vector b, N_Vector x, sunrealtype tol)
Solve the linear system Ax = b to the specificed tolerance.
 Parameters:
b – the righthand side vector
x – the solution vector
tol – the tolerance to solve the system to
 Returns:
gko::LinOp*
the solution

LinearSolver() = default;
10.24. The SUNLINEARSOLVER_KOKKOSDENSE Module
Added in version 6.4.0.
The SUNLINEARSOLVER_KOKKOSDENSE SUNLinearSolver
implementation
provides an interface to KokkosKernels [127] linear solvers
for dense and batched dense (blockdiagonal) systems. Since Kokkos is a modern
C++ library, the module is also written in modern C++ (it requires C++14) as a
header only library. To utilize this SUNLinearSolver
user will need to
include sunlinsol/sunlinsol_kokkosdense.hpp
. More instructions on building
SUNDIALS with Kokkos and KokkosKernels enabled are given in
§1.2.4. For instructions on building and
using Kokkos and KokkosKernels, refer to the
Kokkos
and KokkosKernels.
documentation.
10.24.1. Using SUNLINEARSOLVER_KOKKOSDENSE
The SUNLINEARSOLVER_KOKKOSDENSE module is defined by the DenseLinearSolver
templated class in the sundials::kokkos
namespace:
template<class ExecSpace = Kokkos::DefaultExecutionSpace,
class MemSpace = typename ExecSpace::memory_space>
class DenseLinearSolver : public sundials::impl::BaseLinearSolver,
public sundials::ConvertibleTo<SUNLinearSolver>
To use the SUNLINEARSOLVER_KOKKOSDENSE module, we begin by constructing an instance of a dense linear solver e.g.,
// Create a dense linear solver
sundials::kokkos::DenseLinearSolver<> LS{sunctx};
Instances of the DenseLinearSolver
class are implicitly or explicitly (using
the Convert()
method) convertible to a
SUNLinearSolver
e.g.,
sundials::kokkos::DenseLinearSolver<> LS{sunctx};
SUNLinearSolver LSA = LS; // implicit conversion to SUNLinearSolver
SUNLinearSolver LSB = LS.Convert(); // explicit conversion to SUNLinearSolver
Warning
SUNLinSolFree()
should never be called on a SUNLinearSolver
that
was created via conversion from a sundials::kokkos::DenseLinearSolver
.
Doing so may result in a double free.
The SUNLINEARSOLVER_KOKKOSDENSE module is compatible with the NVECTOR_KOKKOS vector module (see §8.19) and SUNMATRIX_KOKKOSDENSE matrix module (see §9.17).
10.24.2. SUNLINEARSOLVER_KOKKOSDENSE API
In this section we list the public API of the
sundials::kokkos::DenseLinearSolver
class.

template<class ExecSpace = Kokkos::DefaultExecutionSpace, class MemSpace = typename ExecSpace::memory_space>
class DenseLinearSolver : public sundials::impl::BaseLinearSolver, public sundials::ConvertibleTo<SUNLinearSolver> 
DenseLinearSolver() = default;
Default constructor  means the solver must be moved to.

DenseLinearSolver(SUNContext sunctx)
Constructs a new DenseLinearSolver.
 Parameters:
sunctx – The SUNDIALS simulation context (
SUNContext
)

DenseLinearSolver(DenseLinearSolver &&that_solver) noexcept
Move constructor.

DenseLinearSolver &operator=(DenseLinearSolver &&rhs)
Move assignment.

~DenseLinearSolver() override = default
Default destructor.

operator SUNLinearSolver() override
Implicit conversion to a
SUNLinearSolver
.

operator SUNLinearSolver() const override
Implicit conversion to a
SUNLinearSolver
.

SUNLinearSolver Convert() override
Explicit conversion to a
SUNLinearSolver
.

SUNLinearSolver Convert() const override
Explicit conversion to a
SUNLinearSolver
.

DenseLinearSolver() = default;
10.25. SUNLinearSolver Examples
There are SUNLinearSolver
examples that may be installed for each
implementation; these make use of the functions in test_sunlinsol.c
.
These example functions show simple usage of the SUNLinearSolver
family
of modules. The inputs to the examples depend on the linear solver type,
and are output to stdout
if the example is run without the
appropriate number of commandline arguments.
The following is a list of the example functions in test_sunlinsol.c
:
Test_SUNLinSolGetType
: Verifies the returned solver type against the value that should be returned.Test_SUNLinSolGetID
: Verifies the returned solver identifier against the value that should be returned.Test_SUNLinSolInitialize
: Verifies thatSUNLinSolInitialize
can be called and returns successfully.Test_SUNLinSolSetup
: Verifies thatSUNLinSolSetup
can be called and returns successfully.Test_SUNLinSolSolve
: Given aSUNMatrix
object \(A\),N_Vector
objects \(x\) and \(b\) (where \(Ax=b\)) and a desired solution tolerancetol
, this routine clones \(x\) into a new vector \(y\), callsSUNLinSolSolve
to fill \(y\) as the solution to \(Ay=b\) (to the input tolerance), verifies that each entry in \(x\) and \(y\) match to within10*tol
, and overwrites \(x\) with \(y\) prior to returning (in case the calling routine would like to investigate further).Test_SUNLinSolSetATimes
(iterative solvers only): Verifies thatSUNLinSolSetATimes
can be called and returns successfully.Test_SUNLinSolSetPreconditioner
(iterative solvers only): Verifies thatSUNLinSolSetPreconditioner
can be called and returns successfully.Test_SUNLinSolSetScalingVectors
(iterative solvers only): Verifies thatSUNLinSolSetScalingVectors
can be called and returns successfully.Test_SUNLinSolSetZeroGuess
(iterative solvers only): Verifies thatSUNLinSolSetZeroGuess
can be called and returns successfully.Test_SUNLinSolLastFlag
: Verifies thatSUNLinSolLastFlag
can be called, and outputs the result tostdout
.Test_SUNLinSolNumIters
(iterative solvers only): Verifies thatSUNLinSolNumIters
can be called, and outputs the result tostdout
.Test_SUNLinSolResNorm
(iterative solvers only): Verifies thatSUNLinSolResNorm
can be called, and that the result is nonnegative.Test_SUNLinSolResid
(iterative solvers only): Verifies thatSUNLinSolResid
can be called.Test_SUNLinSolSpace
verifies thatSUNLinSolSpace
can be called, and outputs the results tostdout
.
We’ll note that these tests should be performed in a particular
order. For either direct or iterative linear
solvers, Test_SUNLinSolInitialize
must be called
before Test_SUNLinSolSetup
, which must be called
before Test_SUNLinSolSolve
. Additionally, for iterative linear
solvers Test_SUNLinSolSetATimes
, Test_SUNLinSolSetPreconditioner
and Test_SUNLinSolSetScalingVectors
should be called
before Test_SUNLinSolInitialize
;
similarly Test_SUNLinSolNumIters
, Test_SUNLinSolResNorm
and Test_SUNLinSolResid
should be called
after Test_SUNLinSolSolve
. These are called in the appropriate
order in all of the example problems.