11.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 shared-memory N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP or NVECTOR_PTHREADS).

11.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 user-callable 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 §2.4)

Return value:

New SUNLinSol_Band object, or NULL if either A or y are incompatible.

Notes:

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_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.

11.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 in-place 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 bandwidth ml, then the upper triangular factor $U$ can have upper bandwidth as big as smu = MIN(N-1,mu+ml). The lower triangular factor $L$ has lower bandwidth ml.

The SUNLinSol_Band module defines band implementations of all “direct” linear solver operations listed in §11.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 and pivots array to perform the solve.
SUNLinSolLastFlag_Band
SUNLinSolSpace_Band – this only returns information for the storage within the solver object, i.e. storage for N, last_flag, and pivots.
SUNLinSolFree_Band

11.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 shared-memory N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP or NVECTOR_PTHREADS).

11.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 user-callable 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 §2.4)

Return value:

New SUNLinSol_Dense object, or NULL if either A or y are incompatible.

Notes:

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_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.

11.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 in-place 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 §11.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 and pivots array to perform the solve.
SUNLinSolLastFlag_Dense
SUNLinSolSpace_Dense – this only returns information for the storage within the solver object, i.e. storage for N, last_flag, and pivots.
SUNLinSolFree_Dense

11.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 shared-memory N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS).

11.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 user-callable 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 §2.4)

Return value:

New SUNLinSol_KLU object, or NULL if either A or y are incompatible.

Notes:

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_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:
1. 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.
2. 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 previous SUNKLUReInit call).

Return value:

A SUNErrCode

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:
1. AMD,
2. COLAMD, and
3. the natural ordering.
The default is 1 for COLAMD.

Return value:

A SUNErrCode

sun_klu_symbolic *SUNLinSol_KLUGetSymbolic(SUNLinearSolver S)

This function returns a pointer to the KLU symbolic factorization stored in the SUNLinSol_KLU content structure.

When SUNDIALS is compiled with 32-bit indices (SUNDIALS_INDEX_SIZE=32), sun_klu_symbolic is mapped to the KLU type klu_symbolic; when SUNDIALS compiled with 64-bit indices (SUNDIALS_INDEX_SIZE=64) this is mapped to the KLU type klu_l_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.

When SUNDIALS is compiled with 32-bit indices (SUNDIALS_INDEX_SIZE=32), sun_klu_numeric is mapped to the KLU type klu_numeric; when SUNDIALS is compiled with 64-bit indices (SUNDIALS_INDEX_SIZE=64) this is mapped to the KLU type klu_l_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.

When SUNDIALS is compiled with 32-bit indices (SUNDIALS_INDEX_SIZE=32), sun_klu_common is mapped to the KLU type klu_common; when SUNDIALS is compiled with 64-bit indices (SUNDIALS_INDEX_SIZE=64) this is mapped to the KLU type klu_l_common.

11.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 type klu_symbolic or klu_l_symbolic, depending on whether SUNDIALS was installed with 32-bit versus 64-bit indices, respectively,
numeric - KLU storage structure for numeric factorization components, with underlying type klu_numeric or klu_l_numeric, depending on whether SUNDIALS was installed with 32-bit versus 64-bit indices, respectively,
common - storage structure for common KLU solver components, with underlying type klu_common or klu_l_common, depending on whether SUNDIALS was installed with 32-bit versus 64-bit 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 32-bit or 64-bit 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 §2.2.4 for details). Additionally, this wrapper only supports double-precision calculations, and therefore cannot be compiled if SUNDIALS is configured to have sunrealtype set to either extended or single (see §2.3 for details). Since the KLU library supports both 32-bit and 64-bit 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 pre-order the diagonal blocks (the only ones that need to be factored) to reduce fill-in (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 double-precision 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 right-hand side and solution vectors, without requiring costly data copies.

The SUNLinSol_KLU module defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_KLU
SUNLinSolInitialize_KLU – this sets the first_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 integers last_flag and first_factorize. For additional space requirements, see the KLU documentation.
SUNLinSolFree_KLU

11.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 shared-memory N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS). The

11.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 user-callable 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 §2.4)

Return value:

New SUNLinSol_LapackBand object, or NULL if either A or y are incompatible.

Notes:

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_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.

11.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 §2.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 §2.2.4 for details). We note that since there do not exist 128-bit floating-point factorization and solve routines in LAPACK, this interface cannot be compiled when using extended precision for sunrealtype. Similarly, since there do not exist 64-bit 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 in-place 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 bandwidth ml, then the upper triangular factor $U$ can have upper bandwidth as big as smu = MIN(N-1,mu+ml). The lower triangular factor $L$ has lower bandwidth ml.

The SUNLinSol_LapackBand module defines band implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_LapackBand
SUNLinSolInitialize_LapackBand – this does nothing, since all consistency checks are performed at solver creation.
SUNLinSolSetup_LapackBand – this calls either DGBTRF or SGBTRF to perform the $LU$ factorization.
SUNLinSolSolve_LapackBand – this calls either DGBTRS or SGBTRS to use the $LU$ factors and pivots array to perform the solve.
SUNLinSolLastFlag_LapackBand
SUNLinSolSpace_LapackBand – this only returns information for the storage within the solver object, i.e. storage for N, last_flag, and pivots.
SUNLinSolFree_LapackBand

11.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 shared-memory N_Vector implementations (NVECTOR_SERIAL, NVECTOR_OPENMP, or NVECTOR_PTHREADS).

11.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 user-callable 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 §2.4)

Return value:

New SUNLinSol_LapackDense object, or NULL if either A or y are incompatible.

Notes:

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_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.

11.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 §2.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 §2.2.4 for details). We note that since there do not exist 128-bit floating-point factorization and solve routines in LAPACK, this interface cannot be compiled when using extended precision for sunrealtype. Similarly, since there do not exist 64-bit 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 in-place 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 §11.1:

SUNLinSolGetType_LapackDense
SUNLinSolInitialize_LapackDense – this does nothing, since all consistency checks are performed at solver creation.
SUNLinSolSetup_LapackDense – this calls either DGETRF or SGETRF to perform the $LU$ factorization.
SUNLinSolSolve_LapackDense – this calls either DGETRS or SGETRS to use the $LU$ factors and pivots array to perform the solve.
SUNLinSolLastFlag_LapackDense
SUNLinSolSpace_LapackDense – this only returns information for the storage within the solver object, i.e. storage for N, last_flag, and pivots.
SUNLinSolFree_LapackDense

11.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 GPU-enabled 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.

11.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 [122]. The batched LU methods are leveraged when solving block diagonal linear systems of the form

\[\begin{split}\begin{bmatrix} \mathbf{A_0} & 0 & \cdots & 0\\ 0 & \mathbf{A_1} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \mathbf{A_{n-1}}\\ \end{bmatrix} x_j = b_j.\end{split}\]

11.13.2. SUNLinearSolver_MagmaDense Functions

The SUNLinearSolver_MagmaDense module defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_MagmaDense
SUNLinSolInitialize_MagmaDense
SUNLinSolSetup_MagmaDense
SUNLinSolSolve_MagmaDense
SUNLinSolLastFlag_MagmaDense
SUNLinSolFree_MagmaDense

In addition, the module provides the following user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either A or y are incompatible then this routine will return NULL. 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:

A SUNErrCode

11.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;
};

11.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 block-diagonal 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.

11.14.1. SUNLinearSolver_OneMklDense Functions

The SUNLinearSolver_OneMklDense class defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_OneMklDense – returns SUNLINEARSOLVER_ONEMKLDENSE
SUNLinSolInitialize_OneMklDense
SUNLinSolSetup_OneMklDense
SUNLinSolSolve_OneMklDense
SUNLinSolLastFlag_OneMklDense
SUNLinSolFree_OneMklDense

In addition, the class provides the following user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either A or y are incompatible then this routine will return NULL. This routine analyzes the input matrix and vector to determine the linear system size and to assess compatibility with the solver.

11.14.2. SUNLinearSolver_OneMklDense Usage Notes

Warning

The SUNLinearSolver_OneMklDense class only supports 64-bit indexing, thus SUNDIALS must be built for 64-bit 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.

11.15. The SUNLinSol_PCG Module

The SUNLinSol_PCG implementation of the SUNLinearSolver class performs the PCG (Preconditioned Conjugate Gradient [64]) 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$), real-valued 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

(11.8)\[\tilde{A} \tilde{x} = \tilde{b}\]

where

(11.9)\[\begin{split}\tilde{A} &= S P^{-1} A P^{-1} S,\\ \tilde{b} &= S P^{-1} b,\\ \tilde{x} &= S^{-1} P x.\end{split}\]

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:

\[\begin{split}&\| \tilde{b} - \tilde{A} \tilde{x} \|_2 < \delta\\ \Leftrightarrow\quad &\\ &\| S P^{-1} b - S P^{-1} A x \|_2 < \delta\\ \Leftrightarrow\quad &\\ &\| P^{-1} b - P^{-1} A x \|_S < \delta\end{split}\]

where $\| v \|_S = \sqrt{v^T S^T S v}$, with an input tolerance $\delta$.

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either y is incompatible then this routine will return NULL.

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 inputs SUN_PREC_LEFT, SUN_PREC_RIGHT, or SUN_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 one-sided 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:

A SUNErrCode

Notes:

As above, any one of the input values, SUN_PREC_LEFT, SUN_PREC_RIGHT, or SUN_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 non-positive input will result in the default value (5).

Return value:

A SUNErrCode

11.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 most-recent solve,
resnorm - final linear residual norm from the most-recent solve,
last_flag - last error return flag from an internal function,
ATimes - function pointer to perform $Av$ product,
ATData - pointer to structure for ATimes,
Psetup - function pointer to preconditioner setup routine,
Psolve - function pointer to preconditioner solve routine,
PData - pointer to structure for Psetup and Psolve,
s - vector pointer for supplied scaling matrix (default is NULL),
r - a N_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 template N_Vector that is input, and default solver parameters are set.
User-facing “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, and Psolve function pointers and s 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 generic PSetup function and the solver-specific routine (solver-supplied or user-supplied).
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 §11.1:

SUNLinSolGetType_PCG
SUNLinSolInitialize_PCG
SUNLinSolSetATimes_PCG
SUNLinSolSetPreconditioner_PCG
SUNLinSolSetScalingVectors_PCG – since PCG only supports symmetric scaling, the second N_Vector argument to this function is ignored.
SUNLinSolSetZeroGuess_PCG – note the solver assumes a non-zero guess by default and the zero guess flag is reset to SUNFALSE after each call to SUNLinSolSolve_PCG().
SUNLinSolSetup_PCG
SUNLinSolSolve_PCG
SUNLinSolNumIters_PCG
SUNLinSolResNorm_PCG
SUNLinSolResid_PCG
SUNLinSolLastFlag_PCG
SUNLinSolSpace_PCG
SUNLinSolFree_PCG

11.16. The SUNLinSol_SPBCGS Module

The SUNLinSol_SPBCGS implementation of the SUNLinearSolver class performs a Scaled, Preconditioned, Bi-Conjugate Gradient, Stabilized [127] 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.

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either y is incompatible then this routine will return NULL.

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 or SUN_PREC_BOTH the initial guess must be zero (use SUNLinSolSetZeroGuess() 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:

A SUNErrCode

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 non-positive input will result in the default value (5).

Return value:

A SUNErrCode

11.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 most-recent solve,
resnorm - final linear residual norm from the most-recent solve,
last_flag - last error return flag from an internal function,
ATimes - function pointer to perform $Av$ product,
ATData - pointer to structure for ATimes,
Psetup - function pointer to preconditioner setup routine,
Psolve - function pointer to preconditioner solve routine,
PData - pointer to structure for Psetup and Psolve,
s1, s2 - vector pointers for supplied scaling matrices (default is NULL),
r - a N_Vector which holds the current scaled, preconditioned linear system residual,
r_star - a N_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 template N_Vector that is input, and default solver parameters are set.
User-facing “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, and Psolve function pointers and s1 and s2 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 generic PSetup function and the solver-specific routine (solver-supplied or user-supplied).
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 §11.1:

SUNLinSolGetType_SPBCGS
SUNLinSolInitialize_SPBCGS
SUNLinSolSetATimes_SPBCGS
SUNLinSolSetPreconditioner_SPBCGS
SUNLinSolSetScalingVectors_SPBCGS
SUNLinSolSetZeroGuess_SPBCGS – note the solver assumes a non-zero guess by default and the zero guess flag is reset to SUNFALSE after each call to SUNLinSolSolve_SPBCGS().
SUNLinSolSetup_SPBCGS
SUNLinSolSolve_SPBCGS
SUNLinSolNumIters_SPBCGS
SUNLinSolResNorm_SPBCGS
SUNLinSolResid_SPBCGS
SUNLinSolLastFlag_SPBCGS
SUNLinSolSpace_SPBCGS
SUNLinSolFree_SPBCGS

11.17. The SUNLinSol_SPFGMR Module

The SUNLinSol_SPFGMR implementation of the SUNLinearSolver class performs a Scaled, Preconditioned, Flexible, Generalized Minimum Residual [101] 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).

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either y is incompatible then this routine will return NULL.

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 inputs SUN_PREC_LEFT, SUN_PREC_RIGHT, or SUN_PREC_BOTH will result in use of SUN_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 right-preconditioned 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:

A SUNErrCode

Notes:

Since the FGMRES algorithm is designed to only support right preconditioning, then any of the pretype inputs SUN_PREC_LEFT, SUN_PREC_RIGHT, or SUN_PREC_BOTH will result in use of SUN_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 Gram-Schmidt 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:

A SUNErrCode

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:

A SUNErrCode

11.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 Gram-Schmidt orthogonalization (default is modified Gram-Schmidt),
max_restarts - number of FGMRES restarts to allow (default is 0),
numiters - number of iterations from the most-recent solve,
resnorm - final linear residual norm from the most-recent solve,
last_flag - last error return flag from an internal function,
ATimes - function pointer to perform $Av$ product,
ATData - pointer to structure for ATimes,
Psetup - function pointer to preconditioner setup routine,
Psolve - function pointer to preconditioner solve routine,
PData - pointer to structure for Psetup and Psolve,
s1, s2 - vector pointers for supplied scaling matrices (default is NULL),
V - the array of Krylov basis vectors $v_1, \ldots, v_{\text{maxl}+1}$, stored in V[0], ..., V[maxl]. Each $v_i$ is a vector of type N_Vector,
Z - the array of preconditioned Krylov basis vectors $z_1, \ldots, z_{\text{maxl}+1}$, stored in Z[0], ..., Z[maxl]. Each $z_i$ is a vector of type N_Vector,
Hes - the $(\text{maxl}+1)\times\text{maxl}$ Hessenberg matrix. It is stored row-wise so that the (i,j)th element is given by Hes[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 as givens[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 of sunrealtype 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 and vtemp arrays are cloned from a template N_Vector that is input, and default solver parameters are set.
User-facing “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, and Psolve function pointers and s1 and s2 scaling vectors.
In the “initialize” call, the remaining solver data is allocated (V, Hes, givens, and yg )
In the “setup” call, any non-NULL PSetup function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the generic PSetup function and the solver-specific routine (solver-supplied or user-supplied).
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 §11.1:

SUNLinSolGetType_SPFGMR
SUNLinSolInitialize_SPFGMR
SUNLinSolSetATimes_SPFGMR
SUNLinSolSetPreconditioner_SPFGMR
SUNLinSolSetScalingVectors_SPFGMR
SUNLinSolSetZeroGuess_SPFGMR – note the solver assumes a non-zero guess by default and the zero guess flag is reset to SUNFALSE after each call to SUNLinSolSolve_SPFGMR().
SUNLinSolSetup_SPFGMR
SUNLinSolSolve_SPFGMR
SUNLinSolNumIters_SPFGMR
SUNLinSolResNorm_SPFGMR
SUNLinSolResid_SPFGMR
SUNLinSolLastFlag_SPFGMR
SUNLinSolSpace_SPFGMR
SUNLinSolFree_SPFGMR

11.18. The SUNLinSol_SPGMR Module

The SUNLinSol_SPGMR implementation of the SUNLinearSolver class performs a Scaled, Preconditioned, Generalized Minimum Residual [102] 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()).

11.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 user-callable 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 return NULL.

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:

A SUNErrCode

SUNErrCode SUNLinSol_SPGMRSetGSType(SUNLinearSolver S, int gstype)

This function sets the type of Gram-Schmidt 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:

A SUNErrCode

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:

A SUNErrCode

11.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 Gram-Schmidt orthogonalization (default is modified Gram-Schmidt),
max_restarts - number of GMRES restarts to allow (default is 0),
numiters - number of iterations from the most-recent solve,
resnorm - final linear residual norm from the most-recent solve,
last_flag - last error return flag from an internal function,
ATimes - function pointer to perform $Av$ product,
ATData - pointer to structure for ATimes,
Psetup - function pointer to preconditioner setup routine,
Psolve - function pointer to preconditioner solve routine,
PData - pointer to structure for Psetup and Psolve,
s1, s2 - vector pointers for supplied scaling matrices (default is NULL),
V - the array of Krylov basis vectors $v_1, \ldots, v_{\text{maxl}+1}$, stored in V[0], ... V[maxl]. Each $v_i$ is a vector of type N_Vector,
Hes - the $(\text{maxl}+1)\times\text{maxl}$ Hessenberg matrix. It is stored row-wise so that the (i,j)th element is given by Hes[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 as givens[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 of sunrealtype 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 and vtemp arrays are cloned from a template N_Vector that is input, and default solver parameters are set.
User-facing “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, and Psolve function pointers and s1 and s2 scaling vectors.
In the “initialize” call, the remaining solver data is allocated (V, Hes, givens, and yg )
In the “setup” call, any non-NULL PSetup function is called. Typically, this is provided by the SUNDIALS solver itself, that translates between the generic PSetup function and the solver-specific routine (solver-supplied or user-supplied).
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 §11.1:

SUNLinSolGetType_SPGMR
SUNLinSolInitialize_SPGMR
SUNLinSolSetATimes_SPGMR
SUNLinSolSetPreconditioner_SPGMR
SUNLinSolSetScalingVectors_SPGMR
SUNLinSolSetZeroGuess_SPGMR – note the solver assumes a non-zero guess by default and the zero guess flag is reset to SUNFALSE after each call to SUNLinSolSolve_SPGMR().
SUNLinSolSetup_SPGMR
SUNLinSolSolve_SPGMR
SUNLinSolNumIters_SPGMR
SUNLinSolResNorm_SPGMR
SUNLinSolResid_SPGMR
SUNLinSolLastFlag_SPGMR
SUNLinSolSpace_SPGMR
SUNLinSolFree_SPGMR

11.19. The SUNLinSol_SPTFQMR Module

The SUNLinSol_SPTFQMR implementation of the SUNLinearSolver class performs a Scaled, Preconditioned, Transpose-Free Quasi-Minimum Residual [55] 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.

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either y is incompatible then this routine will return NULL.

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 or SUN_PREC_BOTH the initial guess must be zero (use SUNLinSolSetZeroGuess() 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:

A SUNErrCode

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 non-positive input will result in the default value (5).

Return value:

A SUNErrCode

11.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 most-recent solve,
resnorm - final linear residual norm from the most-recent solve,
last_flag - last error return flag from an internal function,
ATimes - function pointer to perform $Av$ product,
ATData - pointer to structure for ATimes,
Psetup - function pointer to preconditioner setup routine,
Psolve - function pointer to preconditioner solve routine,
PData - pointer to structure for Psetup and Psolve,
s1, s2 - vector pointers for supplied scaling matrices (default is NULL),
r_star - a N_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 two N_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 template N_Vector that is input, and default solver parameters are set.
User-facing “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, and Psolve function pointers and s1 and s2 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 generic PSetup function and the solver-specific routine (solver-supplied or user-supplied).
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 §11.1:

SUNLinSolGetType_SPTFQMR
SUNLinSolInitialize_SPTFQMR
SUNLinSolSetATimes_SPTFQMR
SUNLinSolSetPreconditioner_SPTFQMR
SUNLinSolSetScalingVectors_SPTFQMR
SUNLinSolSetZeroGuess_SPTFQMR – note the solver assumes a non-zero guess by default and the zero guess flag is reset to SUNFALSE after each call to SUNLinSolSolve_SPTFQMR().
SUNLinSolSetup_SPTFQMR
SUNLinSolSolve_SPTFQMR
SUNLinSolNumIters_SPTFQMR
SUNLinSolResNorm_SPTFQMR
SUNLinSolResid_SPTFQMR
SUNLinSolLastFlag_SPTFQMR
SUNLinSolSpace_SPTFQMR
SUNLinSolFree_SPTFQMR

11.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).

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object; otherwise this routine will return NULL.

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, and options 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 member Fact 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 is sunrealtype.

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 and U structures. It takes one argument, the SUNLinearSolver 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, the SUNLinearSolver 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.

11.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 grid
lu – pointer to the SuperLU_DIST structure that stores the distributed L and U factors,
scaleperm – pointer to the SuperLU_DIST structure that stores vectors describing the transformations done to the matrix A,
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, 57, 89, 90]. The package uses a SPMD parallel programming model and multithreading to enhance efficiency in distributed-memory 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 §2.2.4 for details). Additionally, the wrapper only supports double-precision 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 32-bit or 64-bit 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 to DOFACT 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 to SamePattern 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 right-hand side and solution vectors, without requiring costly data copies.

The SUNLinSol_SuperLUDIST module defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_SuperLUDIST
SUNLinSolInitialize_SuperLUDIST – this sets the first_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 solves
SUNLinSolSolve_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 integers last_flag and first_factorize. For additional space requirements, see the SuperLU_DIST documentation.
SUNLinSolFree_SuperLUDIST

11.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 shared-memory 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.

11.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 user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object; otherwise this routine will return NULL.

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 and SUNMatrix 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:
1. natural ordering
2. minimal degree ordering on $A^TA$
3. minimal degree ordering on $A^T+A$
4. COLAMD ordering for unsymmetric matrices

The default is 3 for COLAMD.

Return value:

A SUNErrCode

11.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, 88]. 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 §2.2.4 for details). Additionally, this wrapper only supports single- and double-precision calculations, and therefore cannot be compiled if SUNDIALS is configured to have sunrealtype set to extended (see §2.3 for details). Moreover, since the SuperLU_MT library may be installed to support either 32-bit or 64-bit 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 fill-in 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 right-hand side and solution vectors, without requiring costly data copies.

The SUNLinSol_SuperLUMT module defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_SuperLUMT
SUNLinSolInitialize_SuperLUMT – this sets the first_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 integers last_flag and first_factorize. For additional space requirements, see the SuperLU_MT documentation.
SUNLinSolFree_SuperLUMT

11.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.

11.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

\[\begin{split}\begin{bmatrix} \mathbf{A_1} & 0 & \cdots & 0\\ 0 & \mathbf{A_2} & \cdots & 0\\ \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & \cdots & \mathbf{A_n}\\ \end{bmatrix} x_j = b_j\end{split}\]

where all block matrices $\mathbf{A_j}$ share the same sparsity pattern. The matrix must be the SUNMatrix.cuSparse.

11.22.2. SUNLinSol_cuSolverSp_batchQR functions

The SUNLinearSolver_cuSolverSp_batchQR module defines implementations of all “direct” linear solver operations listed in §11.1:

SUNLinSolGetType_cuSolverSp_batchQR
SUNLinSolInitialize_cuSolverSp_batchQR – this sets the first_factorize flag to 1
SUNLinSolSetup_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 system
SUNLinSolSolve_cuSolverSp_batchQR – this calls the cusolverSpXcsrqrsvBatched routine to perform factorization
SUNLinSolLastFlag_cuSolverSp_batchQR
SUNLinSolFree_cuSolverSp_batchQR

In addition, the module provides the following user-callable 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 §2.4)

Return value:

If successful, a SUNLinearSolver object. If either A or y are incompatible then this routine will return NULL.

Notes:

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_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 argument cuSolverInternal and the cuSOLVER batch QR workspace buffer size, in bytes, in the agrument cuSolverWorkspace. 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.

11.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         */
};

11.23. The SUNLINEARSOLVER_GINKGO Module

New 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 §10.16. Instructions on building SUNDIALS with Ginkgo enabled are given in §2.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.

11.23.1. Using SUNLINEARSOLVER_GINKGO

After choosing a compatible N_Vector (see §10.16.1) and creating a Ginkgo-enabled SUNMatrix (see §10.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.

11.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 LinearSolver : public ConvertibleTo<SUNLinearSolver>

LinearSolver() = default;: Default constructor - means the solver must be moved to.

LinearSolver(std::shared_ptr<typename GkoSolverType::Factory> gko_solver_factory, SUNContext sunctx)

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 2-norm.

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 right-hand side vector
x – the solution vector
tol – the tolerance to solve the system to

Returns:

gko::LinOp* the solution

11.24. The SUNLINEARSOLVER_KOKKOSDENSE Module

New in version 6.4.0.

The SUNLINEARSOLVER_KOKKOSDENSE SUNLinearSolver implementation provides an interface to KokkosKernels [123] linear solvers for dense and batched dense (block-diagonal) 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 §2.2.4. For instructions on building and using Kokkos and KokkosKernels, refer to the Kokkos and KokkosKernels. documentation.

11.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 §9.19) and SUNMATRIX_KOKKOSDENSE matrix module (see §10.17).

11.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.

11.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 command-line 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 that SUNLinSolInitialize can be called and returns successfully.
Test_SUNLinSolSetup: Verifies that SUNLinSolSetup can be called and returns successfully.
Test_SUNLinSolSolve: Given a SUNMatrix object $A$, N_Vector objects $x$ and $b$ (where $Ax=b$) and a desired solution tolerance tol, this routine clones $x$ into a new vector $y$, calls SUNLinSolSolve to fill $y$ as the solution to $Ay=b$ (to the input tolerance), verifies that each entry in $x$ and $y$ match to within 10*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 that SUNLinSolSetATimes can be called and returns successfully.
Test_SUNLinSolSetPreconditioner (iterative solvers only): Verifies that SUNLinSolSetPreconditioner can be called and returns successfully.
Test_SUNLinSolSetScalingVectors (iterative solvers only): Verifies that SUNLinSolSetScalingVectors can be called and returns successfully.
Test_SUNLinSolSetZeroGuess (iterative solvers only): Verifies that SUNLinSolSetZeroGuess can be called and returns successfully.
Test_SUNLinSolLastFlag: Verifies that SUNLinSolLastFlag can be called, and outputs the result to stdout.
Test_SUNLinSolNumIters (iterative solvers only): Verifies that SUNLinSolNumIters can be called, and outputs the result to stdout.
Test_SUNLinSolResNorm (iterative solvers only): Verifies that SUNLinSolResNorm can be called, and that the result is non-negative.
Test_SUNLinSolResid (iterative solvers only): Verifies that SUNLinSolResid can be called.
Test_SUNLinSolSpace verifies that SUNLinSolSpace can be called, and outputs the results to stdout.

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.