10.2. ARKODE SUNLinearSolver interface

In Table 10.3, we list the SUNLinSol module functions used within the ARKLS interface. As with the SUNMATRIX module, we emphasize that the ARKODE user does not need to know detailed usage of linear solver functions by the ARKODE code modules in order to use ARKODE. The information is presented as an implementation detail for the interested reader.

Table 10.3 List of SUNLinSol functions called by the ARKODE linear solver interface, depending on the self-identified “type” reported from SUNLinSolGetType(). Functions marked with “X” are required; functions marked with “O” are only called if they are non-NULL in the SUNLinearSolver implementation that is being used.

Routine	DIRECT	ITERATIVE	MATRIX ITERATIVE	MATRIX EMBEDDED
SUNLinSolGetType()	X	X	X	X
SUNLinSolSetATimes()	O	X	O
SUNLinSolSetPreconditioner()	O	O	O
SUNLinSolSetScalingVectors()	O	O	O
SUNLinSolInitialize()	X	X	X
SUNLinSolSetup()	X	X	X
SUNLinSolSolve()	X	X	X	X
SUNLinSolNumIters()1		O	O
SUNLinSolResNorm()2		O	O
SUNLinSolLastFlag()3
SUNLinSolFree()4
SUNLinSolSpace()	O	O	O	O

Notes:

1. SUNLinSolNumIters() is only used to accumulate overall iterative linear solver statistics. If it is not implemented by the SUNLinearSolver module, then ARKLS will consider all solves as requiring zero iterations.

2. Although SUNLinSolResNorm() is optional, if it is not implemented by the SUNLinearSolver then ARKLS will consider all solves a being exact.

3. Although ARKLS does not call SUNLinSolLastFlag() directly, this routine is available for users to query linear solver failure modes.

4. Although ARKLS does not call SUNLinSolFree() directly, this routine should be available for users to call when cleaning up from a simulation.

Since there are a wide range of potential SUNLinSol use cases, the following subsections describe some details of the ARKLS interface, in the case that interested users wish to develop custom SUNLinSol modules.

10.2.1. Lagged matrix information

If the SUNLinSol module identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then it solves a linear system defined by a SUNMATRIX object. ARKLS will update the matrix information infrequently according to the strategies outlined in §2.2.15.2.3. To this end, we differentiate between the desired linear system \(\mathcal A x = b\) with \(\mathcal A = (M-\gamma J)\) and the actual linear system

\[\tilde{\mathcal A} \tilde{x} = b \quad\Leftrightarrow\quad (M-\tilde{\gamma} J)\tilde{x} = b.\]

Since ARKLS updates the SUNMATRIX object infrequently, it is likely that \(\gamma\ne\tilde{\gamma}\), and in turn \(\mathcal A\ne\tilde{\mathcal A}\). Therefore, after calling the SUNLinSol-provided SUNLinSolSolve() routine, we test whether \(\gamma / \tilde{\gamma} \ne 1\), and if this is the case we scale the solution \(\tilde{x}\) to obtain the desired linear system solution \(x\) via

(10.3)\[x = \frac{2}{1 + \gamma / \tilde{\gamma}} \tilde{x}.\]

The motivation for this selection of the scaling factor \(c = 2/(1 + \gamma/\tilde{\gamma})\) follows the derivation in [23, 74]. In short, if we consider a stationary iteration for the linear system as consisting of a solve with \(\tilde{\mathcal A}\) followed with a scaling by \(c\), then for a linear constant-coefficient problem, the error in the solution vector will be reduced at each iteration by the error matrix \(E = I - c \tilde{\mathcal A}^{-1} \mathcal A\), with a convergence rate given by the spectral radius of \(E\). Assuming that stiff systems have a spectrum spread widely over the left half-plane, \(c\) is chosen to minimize the magnitude of the eigenvalues of \(E\).

10.2.2. Iterative linear solver tolerance

If the SUNLinSol object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE, then ARKLS will set the input tolerance delta as described in §2.2.15.3.2. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors() routine is NULL), then ARKLS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

• All solution components have similar magnitude; hence the residual weight vector \(w\) used in the WRMS norm (see §2.2.10), corresponding to the left scaling matrix \(S_1\), should satisfy the assumption

  \[w_i \approx w_{mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

• The SUNLinSol object uses a standard 2-norm to measure convergence.

Under these assumptions, ARKLS adjusts the linear solver convergence requirement as follows (using the notation from (10.2)):

\[\begin{split}&\| \tilde{b} - \tilde{A} \tilde{x} \|_2 < \text{tol}\\ \Leftrightarrow \quad & \| S_1 P_1^{-1} b - S_1 P_1^{-1} A x \|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[w_i \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & w_{mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{w_{mean}}\right)^2\\ \Leftrightarrow \quad & \| P_1^{-1} (b - A x)\|_2 < \frac{\text{tol}}{w_{mean}}\end{split}\]

Therefore we compute the tolerance scaling factor

\[w_{mean} = \|w\|_2 / \sqrt{n}\]

and supply the scaled tolerance delta \(= \text{tol} / w_{mean}\) to the SUNLinSol object.

10.2.3. Providing a custom SUNLinearSolver

In certain instances, users may wish to provide a custom SUNLinSol implementation to ARKODE in order to leverage the structure of a problem. While the “standard” API for these routines is typically sufficient for most users, others may need additional ARKODE-specific information on top of what is provided. For these purposes, we note the following advanced output functions available in ARKStep and MRIStep:

ARKStep advanced outputs: when solving the Newton nonlinear system of equations in predictor-corrector form,

\[\begin{split}\begin{array}{ll} G(z_{cor}) \equiv z_{cor} - \gamma f^I\left(t^I_{n,i}, z_{i} \right) - \tilde{a}_i = 0 &\qquad \text{[$M=I$]},\\ G(z_{cor}) \equiv M z_{cor} - \gamma f^I\left(t^I_{n,i}, z_{i} \right) - \tilde{a}_i = 0 &\qquad \text{[$M$ static]},\\ G(z_{cor}) \equiv M(t^I_{n,i}) (z_{cor} - \tilde{a}_i) - \gamma f^I\left(t^I_{n,i}, z_{i}\right) = 0 &\qquad \text{[$M$ time-dependent]}. \end{array}\end{split}\]

• ARKStepGetCurrentTime() – when called within the computation of a step (i.e., within a solve) this returns \(t^I_{n,i}\). Otherwise the current internal solution time is returned.

• ARKStepGetCurrentState() – when called within the computation of a step (i.e., within a solve) this returns the current stage vector \(z_{i} = z_{cor} + z_{pred}\). Otherwise the current internal solution is returned.

• ARKStepGetCurrentGamma() – returns \(\gamma\).

• ARKStepGetCurrentMassMatrix() – returns \(M(t)\).

• ARKStepGetNonlinearSystemData() – returns \(z_{i}\), \(z_{pred}\), \(f^I(t^I_{n,i}, y_{cur})\), \(\tilde{a}_i\), and \(\gamma\).

MRIStep advanced outputs: when solving the Newton nonlinear system of equations in predictor-corrector form,

\[G(z_{cor}) \equiv z_{cor} - \gamma f^I\left(t^S_{n,i}, z_{i}\right) - \tilde{a}_i = 0\]

• MRIStepGetCurrentTime() – when called within the computation of a step (i.e., within a solve) this returns \(t^S_{n,i}\). Otherwise the current internal solution time is returned.

• MRIStepGetCurrentState() – when called within the computation of a step (i.e., within a solve) this returns the current stage vector \(z_{i} = z_{cor} + z_{pred}\). Otherwise the current internal solution is returned.

• MRIStepGetCurrentGamma() – returns \(\gamma\).

• MRIStepGetNonlinearSystemData() – returns \(z_{i}\), \(z_{pred}\), \(f^I(t^I_{n,i}, y_{cur})\), \(\tilde{a}_i\), and \(\gamma\).

10.3. CVODE SUNLinearSolver interface

Table 10.4 below lists the SUNLinearSolver module linear solver functions used within the CVLS interface. As with the SUNMatrix module, we emphasize that the CVODE user does not need to know detailed usage of linear solver functions by the CVODE code modules in order to use CVODE. The information is presented as an implementation detail for the interested reader.

The linear solver functions listed below are marked with ‘x’ to indicate that they are required, or with \(\dagger\) to indicate that they are only called if they are non-NULL in the SUNLinearSolver implementation that is being used. Note:

1. SUNLinSolNumIters is only used to accumulate overall iterative linear solver statistics. If it is not implemented by the SUNLinearSolver module, then CVLS will consider all solves as requiring zero iterations.

2. Although CVLS does not call SUNLinSolLastFlag directly, this routine is available for users to query linear solver issues directly.

3. Although CVLS does not call SUNLinSolFree directly, this routine should be available for users to call when cleaning up from a simulation.

Table 10.4 List of linear solver function usage in the CVLS interface

	DIRECT	ITERATIVE	MATRIX_ITERATIVE
SUNLinSolGetType()	x	x	x
SUNLinSolSetATimes()	\(\dagger\)	x	\(\dagger\)
SUNLinSolSetPreconditioner()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolSetScalingVectors()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolInitialize()	x	x	x
SUNLinSolSetup()	x	x	x
SUNLinSolSolve()	x	x	x
\(^1\) SUNLinSolNumIters()		\(\dagger\)	\(\dagger\)
\(^2\) SUNLinSolLastFlag()
\(^3\) SUNLinSolFree()
SUNLinSolSpace()	\(\dagger\)	\(\dagger\)	\(\dagger\)

Since there are a wide range of potential SUNLinearSolver use cases, the following subsections describe some details of the CVLS interface, in the case that interested users wish to develop custom SUNLinearSolver modules.

10.3.1. Lagged matrix information

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then the SUNLinearSolver object solves a linear system defined by a SUNMatrix object. CVLS will update the matrix information infrequently according to the strategies outlined in §3.2. To this end, we differentiate between the desired linear system \(Mx=b\) with \(M = (I-\gamma J)\), and the actual linear system

\[\bar{M}\bar{x} = b \quad\Leftrightarrow\quad (I-\bar{\gamma}J)\bar{x} = b.\]

Since CVLS updates the SUNMatrix object infrequently, it is likely that \(\gamma\ne\bar{\gamma}\), and in turn \(M\ne\bar{M}\). When using a BDF method, after calling the SUNLinearSolver-provided SUNLinSolSolve routine, we test whether \(\gamma / \bar{\gamma} \ne 1\), and if this is the case we scale the solution \(\bar{x}\) to correct the linear system solution \(x\) via

(10.4)\[x = \frac{2}{1 + \gamma / \bar{\gamma}} \bar{x}.\]

The motivation for this selection of the scaling factor \(c = 2/(1 +\gamma/\bar{\gamma})\) is discussed in detail in [23, 74]. In short, if we consider a stationary iteration for the linear system as consisting of a solve with \(\bar{M}\) followed by scaling by \(c\), then for a linear constant-coefficient problem, the error in the solution vector will be reduced at each iteration by the error matrix \(E = I - c \bar{M}^{-1} M\), with a convergence rate given by the spectral radius of \(E\). Assuming that stiff systems have a spectrum spread widely over the left half-plane, \(c\) is chosen to minimize the magnitude of the eigenvalues of \(E\).

10.3.2. Iterative linear solver tolerance

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE then CVLS will set the input tolerance delta as described in §3.2.1. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors routine is NULL), then CVLS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

1. All solution components have similar magnitude; hence the error weight vector \(W\) used in the WRMS norm (see §3.2.1) should satisfy the assumption

   \[W_i \approx W_{mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

2. The SUNLinearSolver object uses a standard 2-norm to measure convergence.

Since CVODE uses identical left and right scaling matrices, \(S_1 = S_2 = S = \operatorname{diag}(W)\), then the linear solver convergence requirement is converted as follows (using the notation from equations (10.1) – (10.2)):

\[\begin{split}\begin{aligned} &\left\| \tilde{b} - \tilde{A} \tilde{x} \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \left\| S P_1^{-1} b - S P_1^{-1} A x \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[W_i \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & W_{mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{W_{mean}}\right)^2\\ \Leftrightarrow \quad & \left\| P_1^{-1} (b - A x)\right\|_2 < \frac{\text{tol}}{W_{mean}}\end{aligned}\end{split}\]

Therefore the tolerance scaling factor

\[W_{mean} = \|W\|_2 / \sqrt{n}\]

is computed and the scaled tolerance delta\(= \text{tol} / W_{mean}\) is supplied to the SUNLinearSolver object.

10.4. CVODES SUNLinearSolver interface

Table 10.5 below lists the SUNLinearSolver module linear solver functions used within the CVLS interface. As with the SUNMatrix module, we emphasize that the CVODES user does not need to know detailed usage of linear solver functions by the CVODES code modules in order to use CVODES. The information is presented as an implementation detail for the interested reader.

The linear solver functions listed below are marked with “x” to indicate that they are required, or with “\(\dagger\)” to indicate that they are only called if they are non-NULL in the SUNLinearSolver implementation that is being used. Note:

1. SUNLinSolNumIters is only used to accumulate overall iterative linear solver statistics. If it is not implemented by the SUNLinearSolver module, then CVLS will consider all solves as requiring zero iterations.

2. Although CVLS does not call SUNLinSolLastFlag directly, this routine is available for users to query linear solver issues directly.

3. Although CVLS does not call SUNLinSolFree directly, this routine should be available for users to call when cleaning up from a simulation.

Table 10.5 List of linear solver function usage in the CVLS interface

	DIRECT	ITERATIVE	MATRIX_ITERATIVE
SUNLinSolGetType()	x	x	x
SUNLinSolSetATimes()	\(\dagger\)	x	\(\dagger\)
SUNLinSolSetPreconditioner()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolSetScalingVectors()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolInitialize()	x	x	x
SUNLinSolSetup()	x	x	x
SUNLinSolSolve()	x	x	x
\(^1\) SUNLinSolNumIters()		\(\dagger\)	\(\dagger\)
\(^2\) SUNLinSolLastFlag()
\(^3\) SUNLinSolFree()
SUNLinSolSpace()	\(\dagger\)	\(\dagger\)	\(\dagger\)

Since there are a wide range of potential SUNLinearSolver use cases, the following subsections describe some details of the CVLS interface, in the case that interested users wish to develop custom SUNLinearSolver modules.

10.4.1. Lagged matrix information

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then the SUNLinearSolver object solves a linear system defined by a SUNMatrix object. CVLS will update the matrix information infrequently according to the strategies outlined in §4.2. To this end, we differentiate between the desired linear system \(Mx=b\) with \(M = (I-\gamma J)\), and the actual linear system

\[\bar{M}\bar{x} = b \quad\Leftrightarrow\quad (I-\bar{\gamma}J)\bar{x} = b.\]

Since CVLS updates the SUNMatrix object infrequently, it is likely that \(\gamma\ne\bar{\gamma}\), and in turn \(M\ne\bar{M}\). When using a BDF method, after calling the SUNLinearSolver-provided SUNLinSolSolve routine, we test whether \(\gamma / \bar{\gamma} \ne 1\), and if this is the case we scale the solution \(\bar{x}\) to correct the linear system solution \(x\) via

(10.5)\[x = \frac{2}{1 + \gamma / \bar{\gamma}} \bar{x}.\]

The motivation for this selection of the scaling factor \(c = 2/(1 +\gamma/\bar{\gamma})\) is discussed in detail in [23, 74]. In short, if we consider a stationary iteration for the linear system as consisting of a solve with \(\bar{M}\) followed by scaling by \(c\), then for a linear constant-coefficient problem, the error in the solution vector will be reduced at each iteration by the error matrix \(E = I - c \bar{M}^{-1} M\), with a convergence rate given by the spectral radius of \(E\). Assuming that stiff systems have a spectrum spread widely over the left half-plane, \(c\) is chosen to minimize the magnitude of the eigenvalues of \(E\).

10.4.2. Iterative linear solver tolerance

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE then CVLS will set the input tolerance delta as described in §4.2.1. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors routine is NULL), then CVLS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

1. All solution components have similar magnitude; hence the error weight vector \(W\) used in the WRMS norm (see §4.2.1) should satisfy the assumption

   \[W_i \approx W_{mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

2. The SUNLinearSolver object uses a standard 2-norm to measure convergence.

Since CVODES uses identical left and right scaling matrices, \(S_1 = S_2 = S = \operatorname{diag}(W)\), then the linear solver convergence requirement is converted as follows (using the notation from equations (10.1) – (10.2)):

\[\begin{split}\begin{aligned} &\| \tilde{b} - \tilde{A} \tilde{x} \|_2 < \text{tol}\\ \Leftrightarrow \quad & \| S P_1^{-1} b - S P_1^{-1} A x \|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[W_i \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & W_{mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{W_{mean}}\right)^2\\ \Leftrightarrow \quad & \| P_1^{-1} (b - A x)\|_2 < \frac{\text{tol}}{W_{mean}}\end{aligned}\end{split}\]

Therefore the tolerance scaling factor

\[W_{mean} = \|W\|_2 / \sqrt{n}\]

is computed and the scaled tolerance delta\(= \text{tol} / W_{mean}\) is supplied to the SUNLinearSolver object.

10.5. IDA SUNLinearSolver interface

Table 10.6 below lists the SUNLinearSolver module linear solver functions used within the IDALS interface. As with the SUNMatrix module, we emphasize that the IDA user does not need to know detailed usage of linear solver functions by the IDA code modules in order to use IDA. The information is presented as an implementation detail for the interested reader.

The linear solver functions listed below are marked with ‘x’ to indicate that they are required, or with \(\dagger\) to indicate that they are only called if they are non-NULL in the SUNLinearSolver implementation that is being used. Note:

1. Although IDALS does not call SUNLinSolLastFlag directly, this routine is available for users to query linear solver issues directly.

2. Although IDALS does not call SUNLinSolFree directly, this routine should be available for users to call when cleaning up from a simulation.

Table 10.6 List of linear solver function usage in the IDALS interface

	DIRECT	ITERATIVE	MATRIX_ITERATIVE
SUNLinSolGetType()	x	x	x
SUNLinSolSetATimes()	\(\dagger\)	x	\(\dagger\)
SUNLinSolSetPreconditioner()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolSetScalingVectors()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolInitialize()	x	x	x
SUNLinSolSetup()	x	x	x
SUNLinSolSolve()	x	x	x
SUNLinSolNumIters()		x	x
SUNLinSolResid()		x	x
\(^1\) SUNLinSolLastFlag()
\(^2\)SUNLinSolFree()
SUNLinSolSpace()	\(\dagger\)	\(\dagger\)	\(\dagger\)

Since there are a wide range of potential SUNLinearSolver use cases, the following subsections describe some details of the IDALS interface, in the case that interested users wish to develop custom SUNLinearSolver modules.

10.5.1. Lagged matrix information

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then the SUNLinearSolver object solves a linear system defined by a SUNMatrix object. IDALS will update the matrix information infrequently according to the strategies outlined in §5.2. To this end, we differentiate between the desired linear system \(Jx=b\) with \(J = \left(\dfrac{\partial F}{\partial y}-c_j \dfrac{\partial F}{\partial\dot{y}}\right)\), and the actual linear system \(\bar{J}\bar{x}=b\) with

\[\bar{J} = \dfrac{\partial \bar{F}}{\partial y} - \bar{c}_j \dfrac{\partial \bar{F}}{\partial\dot{y}},\]

where the overlines indicate the lagged versions of these numbers and matrices.

Since IDALS updates the SUNMatrix objects infrequently and it is likely that \(c_j\ne\bar{c}_j\), then typically \(J\ne\bar{J}\). Thus after calling the SUNLinearSolver-provided SUNLinSolSolve routine, we test whether \(\dfrac{c_j}{\bar{c}_j} \ne 1\), and if this is the case we scale the solution \(\bar{x}\) to correct the linear system solution \(x\) via

(10.6)\[x = \frac{2}{1 + c_j / \bar{c}_j} \bar{x}.\]

The motivation for this selection of the scaling factor \(c = 2/(1 + c_j/\bar{c}_j)\) is discussed in detail in [23, 74]. In short, if we consider a stationary iteration for the linear system as consisting of a solve with \(\bar{J}\) followed by scaling by \(c\), then for a linear constant-coefficient problem, the error in the solution vector will be reduced at each iteration by the error matrix \(E = I - c \bar{J}^{-1} J\), with a convergence rate given by the spectral radius of \(E\). Assuming that stiff systems have a spectrum spread widely over the left half-plane, \(c\) is chosen to minimize the magnitude of the eigenvalues of \(E\).

10.5.2. Iterative linear solver tolerance

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE then IDALS will set the input tolerance delta as described in §5.2.2. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors routine is NULL), then IDALS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

1. All solution components have similar magnitude; hence the error weight vector \(W\) used in the WRMS norm (see §5.2.2) should satisfy the assumption

   \[W_i \approx W_{mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

2. The SUNLinearSolver object uses a standard 2-norm to measure convergence.

Since IDA uses identical left and right scaling matrices, \(S_1 = S_2 = S = \operatorname{diag}(W)\), then the linear solver convergence requirement is converted as follows (using the notation from equations (10.1) – (10.2)):

\[\begin{split}\begin{aligned} &\left\| \tilde{b} - \tilde{A} \tilde{x} \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \left\| S P_1^{-1} b - S P_1^{-1} A x \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[W_i \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & W_{mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{W_{mean}}\right)^2\\ \Leftrightarrow \quad & \left\| P_1^{-1} (b - A x)\right\|_2 < \frac{\text{tol}}{W_{mean}}\end{aligned}\end{split}\]

Therefore the tolerance scaling factor

\[W_{mean} = \|W\|_2 / \sqrt{n}\]

is computed and the scaled tolerance delta\(= \text{tol} / W_{mean}\) is supplied to the SUNLinearSolver object.

10.6. IDAS SUNLinearSolver interface

Table 10.7 below lists the SUNLinearSolver module linear solver functions used within the IDALS interface. As with the SUNMatrix module, we emphasize that the IDA user does not need to know detailed usage of linear solver functions by the IDA code modules in order to use IDA. The information is presented as an implementation detail for the interested reader.

The linear solver functions listed below are marked with ‘x’ to indicate that they are required, or with \(\dagger\) to indicate that they are only called if they are non-NULL in the SUNLinearSolver implementation that is being used. Note:

1. Although IDALS does not call SUNLinSolLastFlag directly, this routine is available for users to query linear solver issues directly.

2. Although IDALS does not call SUNLinSolFree directly, this routine should be available for users to call when cleaning up from a simulation.

Table 10.7 List of linear solver function usage in the IDALS interface

	DIRECT	ITERATIVE	MATRIX_ITERATIVE
SUNLinSolGetType()	x	x	x
SUNLinSolSetATimes()	\(\dagger\)	x	\(\dagger\)
SUNLinSolSetPreconditioner()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolSetScalingVectors()	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolInitialize()	x	x	x
SUNLinSolSetup()	x	x	x
SUNLinSolSolve()	x	x	x
SUNLinSolNumIters()		x	x
SUNLinSolResid()		x	x
\(^1\) SUNLinSolLastFlag()
\(^2\)SUNLinSolFree()
SUNLinSolSpace()	\(\dagger\)	\(\dagger\)	\(\dagger\)

Since there are a wide range of potential SUNLinearSolver use cases, the following subsections describe some details of the IDALS interface, in the case that interested users wish to develop custom SUNLinearSolver modules.

10.6.1. Lagged matrix information

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then the SUNLinearSolver object solves a linear system defined by a SUNMatrix object. IDALS will update the matrix information infrequently according to the strategies outlined in §6.2. To this end, we differentiate between the desired linear system \(Jx=b\) with \(J = \left(\dfrac{\partial F}{\partial y}-c_j \dfrac{\partial F}{\partial\dot{y}}\right)\), and the actual linear system \(\bar{J}\bar{x}=b\) with

\[\bar{J} = \dfrac{\partial \bar{F}}{\partial y} - \bar{c}_j \dfrac{\partial \bar{F}}{\partial\dot{y}},\]

where the overlines indicate the lagged versions of these numbers and matrices.

Since IDALS updates the SUNMatrix objects infrequently and it is likely that \(c_j\ne\bar{c}_j\), then typically \(J\ne\bar{J}\). Thus after calling the SUNLinearSolver-provided SUNLinSolSolve routine, we test whether \(\dfrac{c_j}{\bar{c}_j} \ne 1\), and if this is the case we scale the solution \(\bar{x}\) to correct the linear system solution \(x\) via

(10.7)\[x = \frac{2}{1 + c_j / \bar{c}_j} \bar{x}.\]

The motivation for this selection of the scaling factor \(c = 2/(1 + c_j/\bar{c}_j)\) is discussed in detail in [23, 74]. In short, if we consider a stationary iteration for the linear system as consisting of a solve with \(\bar{J}\) followed by scaling by \(c\), then for a linear constant-coefficient problem, the error in the solution vector will be reduced at each iteration by the error matrix \(E = I - c \bar{J}^{-1} J\), with a convergence rate given by the spectral radius of \(E\). Assuming that stiff systems have a spectrum spread widely over the left half-plane, \(c\) is chosen to minimize the magnitude of the eigenvalues of \(E\).

10.6.2. Iterative linear solver tolerance

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE then IDALS will set the input tolerance delta as described in §6.2.2. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors routine is NULL), then IDALS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

1. All solution components have similar magnitude; hence the error weight vector \(W\) used in the WRMS norm (see §6.2.2) should satisfy the assumption

   \[W_i \approx W_{mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

2. The SUNLinearSolver object uses a standard 2-norm to measure convergence.

Since IDA uses identical left and right scaling matrices, \(S_1 = S_2 = S = \operatorname{diag}(W)\), then the linear solver convergence requirement is converted as follows (using the notation from equations (10.1) – (10.2)):

\[\begin{split}\begin{aligned} &\left\| \tilde{b} - \tilde{A} \tilde{x} \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \left\| S P_1^{-1} b - S P_1^{-1} A x \right\|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[W_i \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & W_{mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{W_{mean}}\right)^2\\ \Leftrightarrow \quad & \left\| P_1^{-1} (b - A x)\right\|_2 < \frac{\text{tol}}{W_{mean}}\end{aligned}\end{split}\]

Therefore the tolerance scaling factor

\[W_{mean} = \|W\|_2 / \sqrt{n}\]

is computed and the scaled tolerance delta\(= \text{tol} / W_{mean}\) is supplied to the SUNLinearSolver object.

10.7. KINSOL SUNLinearSolver interface

Table 10.8 below lists the SUNLinearSolver module linear solver functions used within the KINLS interface. As with the SUNMatrix module, we emphasize that the KINSOL user does not need to know detailed usage of linear solver functions by the KINSOL code modules in order to use KINSOL. The information is presented as an implementation detail for the interested reader.

The linear solver functions listed below are marked with “x” to indicate that they are required, or with “\(\dagger\)” to indicate that they are only called if they are non-NULL in the SUNLinearSolver implementation that is being used. Note:

1. SUNLinSolNumIters() is only used to accumulate overall iterative linear solver statistics. If it is not implemented by the SUNLinearSolver module, then KINLS will consider all solves as requiring zero iterations.

2. Although SUNLinSolResNorm() is optional, if it is not implemented by the SUNLinearSolver then KINLS will consider all solves a being exact.

3. Although KINLS does not call SUNLinSolLastFlag() directly, this routine is available for users to query linear solver issues directly.

4. Although KINLS does not call SUNLinSolFree() directly, this routine should be available for users to call when cleaning up from a simulation.

Table 10.8 List of linear solver function usage in the KINLS interface

SUNLinSolGetType	x	x	x
SUNLinSolSetATimes	\(\dagger\)	x	\(\dagger\)
SUNLinSolSetPreconditioner	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolSetScalingVectors	\(\dagger\)	\(\dagger\)	\(\dagger\)
SUNLinSolInitialize	x	x	x
SUNLinSolSetup	x	x	x
SUNLinSolSolve	x	x	x
\(^1\)SUNLinSolNumIters		\(\dagger\)	\(\dagger\)
\(^2\)SUNLinSolResNorm		\(\dagger\)	\(\dagger\)
\(^3\)SUNLinSolLastFlag
\(^4\)SUNLinSolFree
SUNLinSolSpace	\(\dagger\)	\(\dagger\)	\(\dagger\)

Since there are a wide range of potential SUNLinearSolver use cases, the following subsections describe some details of the KINLS interface, in the case that interested users wish to develop custom SUNLinearSolver modules.

10.7.1. Lagged matrix information

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_DIRECT or SUNLINEARSOLVER_MATRIX_ITERATIVE, then the SUNLinearSolver object solves a linear system defined by a SUNMatrix object. As a result, KINSOL can perform its optional residual monitoring scheme, described in §7.2.8.

10.7.2. Iterative linear solver tolerance

If the SUNLinearSolver object self-identifies as having type SUNLINEARSOLVER_ITERATIVE or SUNLINEARSOLVER_MATRIX_ITERATIVE then KINLS will adjust the linear solver tolerance delta as described in §7.2.9 during the course of the nonlinear solve process. However, if the iterative linear solver does not support scaling matrices (i.e., the SUNLinSolSetScalingVectors routine is NULL), then KINLS will be unable to fully handle ill-conditioning in the nonlinear solve process through the solution and residual scaling operators described in §7.2.4. In this case, KINLS will attempt to adjust the linear solver tolerance to account for this lack of functionality. To this end, the following assumptions are made:

1. All residual components have similar magnitude; hence the scaling matrix \(D_F\) used in computing the linear residual norm (see §7.2.4) should satisfy the assumption

   \[(D_F)_{i,i} \approx D_{F,mean},\quad \text{for}\quad i=0,\ldots,n-1.\]

2. The SUNLinearSolver object uses a standard 2-norm to measure convergence.

Since KINSOL uses \(D_F\) as the left-scaling matrix, \(S_1 = D_F\), then the linear solver convergence requirement is converted as follows (using the notation from equations (10.1) – (10.2):

\[\begin{split}&\| \tilde{b} - \tilde{A} \tilde{x} \|_2 < \text{tol}\\ \Leftrightarrow \quad & \| D_F P_1^{-1} b - D_F P_1^{-1} A x \|_2 < \text{tol}\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[(D_F)_{i,i} \left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & D_{F,mean}^2 \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \text{tol}^2\\ \Leftrightarrow \quad & \sum_{i=0}^{n-1} \left[\left(P_1^{-1} (b - A x)\right)_i\right]^2 < \left(\frac{\text{tol}}{D_{F,mean}}\right)^2\\ \Leftrightarrow \quad & \| P_1^{-1} (b - A x)\|_2 < \frac{\text{tol}}{D_{F,mean}}\end{split}\]

Therefore the tolerance scaling factor

\[D_{F,mean} = \frac{1}{\sqrt{n}}\left(\sum_{i=0}^{n-1} (D_F)_{i,i}^2\right)^{1/2}\]

is computed and the scaled tolerance delta\(= \text{tol} / D_{F,mean}\) is supplied to the SUNLinearSolver object.

10.7.3. Matrix-embedded solver incompatibility

At present, KINLS is incompatible with SUNLinearSolver objects that self-identify as having type SUNLINEARSOLVER_MATRIX_EMBEDDED. Support for such user-supplied linear solvers may be added in a future release. Users interested in such support are recommended to contact the SUNDIALS development team.