2.4.4. Usersupplied functions
The usersupplied functions for ARKODE consist of:
at least one function defining the ODE (required),
a function that provides the error weight vector (optional),
a function that provides the residual weight vector (optional),
a function that handles explicit time step stability (optional),
a function that updates the implicit stage prediction (optional),
a function that defines auxiliary temporal rootfinding problem(s) to solve (optional),
one or two functions that provide Jacobianrelated information for the linear solver, if a component is treated implicitly and a Newtonbased nonlinear iteration is chosen (optional),
one or two functions that define the preconditioner for use in any of the Krylov iterative algorithms, if linear systems of equations are to be solved using an iterative method (optional),
if the problem involves a nonidentity mass matrix \(M\ne I\) with ARKStep:
one or two functions that provide massmatrixrelated information for the linear and mass matrix solvers (required),
one or two functions that define the mass matrix preconditioner for use if an iterative mass matrix solver is chosen (optional), and
a function that handles vector resizing operations, if the underlying vector structure supports resizing (as opposed to deletion/recreation), and if the user plans to call
ARKodeResize()
(optional).MRIStep only: functions to be called before and after each inner integration to perform any communication or memory transfers of forcing data supplied by the outer integrator to the inner integrator, or state data supplied by the inner integrator to the outer integrator.
if relaxation is enabled (optional), a function that evaluates the conservative or dissipative function \(\xi(y(t))\) (required) and a function to evaluate its Jacobian \(\xi'(y(t))\) (required).
2.4.4.1. ODE righthand side
The user must supply at least one function of type ARKRhsFn
to
specify the IVPdefininig righthand side function(s) when creating the
ARKODE timestepping module:

typedef int (*ARKRhsFn)(sunrealtype t, N_Vector y, N_Vector ydot, void *user_data)
These functions compute the ODE righthand side for a given value of the independent variable \(t\) and state vector \(y\).
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param ydot:
the output vector that forms [a portion of] the ODE RHS \(f(t,y)\).
 Param user_data:
the user_data pointer that was passed to
ARKodeSetUserData()
. Return:
An ARKRhsFn should return 0 if successful, a positive value if a recoverable error occurred (in which case ARKODE will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and ARK_RHSFUNC_FAIL is returned).
Note
Allocation of memory for ydot is handled within ARKODE.
The vector ydot may be uninitialized on input; it is the user’s responsibility to fill this entire vector with meaningful values.
A recoverable failure error return from the ARKRhsFn is typically used to flag a value of the dependent variable \(y\) that is “illegal” in some way (e.g., negative where only a nonnegative value is physically meaningful). If such a return is made, ARKODE will attempt to recover (possibly repeating the nonlinear iteration, or reducing the step size in ARKodeEvolve) in order to avoid this recoverable error return. There are some situations in which recovery is not possible even if the righthand side function returns a recoverable error flag. One is when this occurs at the very first call to the ARKRhsFn (in which case ARKODE returns ARK_FIRST_RHSFUNC_ERR). Another is when a recoverable error is reported by ARKRhsFn after the timestepping module completes a successful stage, in which case ARKodeEvolve returns ARK_UNREC_RHSFUNC_ERR). Finally, when ARKODE is run in fixedstep mode, it may halt on a recoverable error flag that would normally have resulted in a stepsize reduction.
2.4.4.2. Error weight function
As an alternative to providing the relative and absolute tolerances,
the user may provide a function of type ARKEwtFn
to compute a
vector ewt containing the weights in the WRMS norm
\(\v\_{WRMS} = \left(\dfrac{1}{n} \displaystyle \sum_{i=1}^n \left(ewt_i\; v_i\right)^2
\right)^{1/2}\). These weights will be used in place of those defined
in §2.2.7.

typedef int (*ARKEwtFn)(N_Vector y, N_Vector ewt, void *user_data)
This function computes the WRMS error weights for the vector \(y\).
 Param y:
the dependent variable vector at which the weight vector is to be computed.
 Param ewt:
the output vector containing the error weights.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
function Return:
An ARKEwtFn function must return 0 if it successfully set the error weights, and 1 otherwise.
Note
Allocation of memory for ewt is handled within ARKODE.
The error weight vector must have all components positive. It is the user’s responsibility to perform this test and return 1 if it is not satisfied.
2.4.4.3. Residual weight function
Warning
The functions in this section are specific to timestepping modules that support nonidentity mass matrices.
As an alternative to providing the scalar or vector absolute residual
tolerances (when the IVP units differ from the solution units), the
user may provide a function of type ARKRwtFn
to compute a
vector rwt containing the weights in the WRMS norm
\(\v\_{WRMS} = \left(\dfrac{1}{n} \displaystyle \sum_{i=1}^n \left(rwt_i\; v_i\right)^2
\right)^{1/2}\). These weights will be used in place of those defined
in §2.2.7.

typedef int (*ARKRwtFn)(N_Vector y, N_Vector rwt, void *user_data)
This function computes the WRMS residual weights for the vector \(y\).
 Param y:
the dependent variable vector at which the weight vector is to be computed.
 Param rwt:
the output vector containing the residual weights.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
An ARKRwtFn function must return 0 if it successfully set the residual weights, and 1 otherwise.
Note
Allocation of memory for rwt is handled within ARKODE.
The residual weight vector must have all components positive. It is the user’s responsibility to perform this test and return 1 if it is not satisfied.
2.4.4.4. Time step adaptivity function
Warning
The function in this section is only used in nowdeprecated functions in ARKStep and ERKStep, and will be removed in a future release.
As an alternative to using one of the builtin time step adaptivity
methods for controlling solution error, the user may provide a
function of type ARKAdaptFn
to compute a target step size
\(h\) for the next integration step. These steps should be chosen
such that the error estimate for the next time step remains below 1.

typedef int (*ARKAdaptFn)(N_Vector y, sunrealtype t, sunrealtype h1, sunrealtype h2, sunrealtype h3, sunrealtype e1, sunrealtype e2, sunrealtype e3, int q, int p, sunrealtype *hnew, void *user_data)
This function implements a time step adaptivity algorithm that chooses \(h\) to satisfy the error tolerances.
 Param y:
the current value of the dependent variable vector.
 Param t:
the current value of the independent variable.
 Param h1:
the current step size, \(t_n  t_{n1}\).
 Param h2:
the previous step size, \(t_{n1}  t_{n2}\).
 Param h3:
the step size \(t_{n2}t_{n3}\).
 Param e1:
the error estimate from the current step, \(n\).
 Param e2:
the error estimate from the previous step, \(n1\).
 Param e3:
the error estimate from the step \(n2\).
 Param q:
the global order of accuracy for the method.
 Param p:
the global order of accuracy for the embedded method.
 Param hnew:
the output value of the next step size.
 Param user_data:
a pointer to user data, the same as the h_data parameter that was passed to
ARKStepSetAdaptivityFn()
orERKStepSetAdaptivityFn()
. Return:
An ARKAdaptFn function should return 0 if it successfully set the next step size, and a nonzero value otherwise.
Deprecated since version 5.7.0: Use the SUNAdaptController infrastructure instead (see §12.1).
2.4.4.5. Explicit stability function
Warning
The functions in this section are specific to timestepping modules that support temporal adaptivity.
A user may supply a function to predict the maximum stable step size
for an explicit portion of their IVP. While the accuracybased time step
adaptivity algorithms may be sufficient for retaining a stable solution to
the ODE system, these may be inefficient if the explicit righthand side
function contains moderately stiff terms. In this scenario, a user may
provide a function of type ARKExpStabFn
to provide this stability information to ARKODE. This function
must set the scalar step size satisfying the stability restriction for
the upcoming time step. This value will subsequently be bounded by
the usersupplied values for the minimum and maximum allowed time
step, and the accuracybased time step.

typedef int (*ARKExpStabFn)(N_Vector y, sunrealtype t, sunrealtype *hstab, void *user_data)
This function predicts the maximum stable step size for the explicit portion of the ODE system.
 Param y:
the current value of the dependent variable vector.
 Param t:
the current value of the independent variable.
 Param hstab:
the output value with the absolute value of the maximum stable step size.
 Param user_data:
a pointer to user data, the same as the estab_data parameter that was passed to
ARKodeSetStabilityFn()
. Return:
An ARKExpStabFn function should return 0 if it successfully set the upcoming stable step size, and a nonzero value otherwise.
Note
If this function is not supplied, or if it returns hstab \(\le 0.0\), then ARKODE will assume that there is no explicit stability restriction on the time step size.
2.4.4.6. Implicit stage prediction function
A user may supply a function to update the prediction for each implicit stage solution.
If supplied, this routine will be called after any existing ARKODE predictor
algorithm completes, so that the predictor may be modified by the user as desired.
In this scenario, a user may provide a function of type ARKStagePredictFn
to provide this implicit predictor to ARKODE. This function takes as input the
alreadypredicted implicit stage solution and the corresponding “time” for that prediction;
it then updates the prediction vector as desired. If the usersupplied routine will
construct a full prediction (and thus the ARKODE prediction is irrelevant), it is
recommended that the user not call ARKodeSetPredictorMethod()
, thereby leaving
the default trivial predictor in place.

typedef int (*ARKStagePredictFn)(sunrealtype t, N_Vector zpred, void *user_data)
This function updates the prediction for the implicit stage solution.
 Param t:
the current value of the independent variable containing the “time” corresponding to the predicted solution.
 Param zpred:
the ARKODEpredicted stage solution on input, and the usermodified predicted stage solution on output.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
An ARKStagePredictFn function should return 0 if it successfully set the upcoming stable step size, and a nonzero value otherwise.
Note
This may be useful if there are bound constraints on the solution, and these should be enforced prior to beginning the nonlinear or linear implicit solver algorithm.
This routine is incompatible with the “minimum correction predictor” – option 5 to the routine
ARKodeSetPredictorMethod()
. If both are selected, then ARKODE will override its builtin implicit predictor routine to instead use option 0 (trivial predictor).
2.4.4.7. Rootfinding function
If a rootfinding problem is to be solved during integration of the
ODE system, the user must supply a function of type ARKRootFn
.

typedef int (*ARKRootFn)(sunrealtype t, N_Vector y, sunrealtype *gout, void *user_data)
This function implements a vectorvalued function \(g(t,y)\) such that roots are sought for the components \(g_i(t,y)\), \(i=0,\ldots,\) nrtfn1.
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param gout:
the output array, of length nrtfn, with components \(g_i(t,y)\).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to the
SetUserData
function Return:
An ARKRootFn function should return 0 if successful or a nonzero value if an error occurred (in which case the integration is halted and ARKODE returns ARK_RTFUNC_FAIL).
Note
Allocation of memory for gout is handled within ARKODE.
2.4.4.8. Jacobian construction
If a matrixbased linear solver module is used (i.e., a nonNULL SUNMatrix
object was supplied to ARKodeSetLinearSolver()
, the user may provide a
function of type ARKLsJacFn
to provide the Jacobian approximation or
ARKLsLinSysFn
to provide an approximation of the linear system
\(\mathcal{A}(t,y) = M(t)  \gamma J(t,y)\).

typedef int (*ARKLsJacFn)(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
This function computes the Jacobian matrix \(J(t,y) = \dfrac{\partial f^I}{\partial y}(t,y)\) (or an approximation to it).
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector, namely the predicted value of \(y(t)\).
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param Jac:
the output Jacobian matrix.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Param tmp*:
pointers to memory allocated to variables of type
N_Vector
which can be used by an ARKLsJacFn as temporary storage or work space. Return:
An ARKLsJacFn function should return 0 if successful, a positive value if a recoverable error occurred (in which case ARKODE will attempt to correct, while ARKLS sets last_flag to ARKLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted,
ARKodeEvolve()
returns ARK_LSETUP_FAIL and ARKLS sets last_flag to ARKLS_JACFUNC_UNRECVR).
Note
Information regarding the specific
SUNMatrix
structure (e.g.~number of rows, upper/lower bandwidth, sparsity type) may be obtained through using the implementationspecificSUNMatrix
interface functions (see §9 for details).When using a linear solver of type
SUNLINEARSOLVER_DIRECT
, prior to calling the usersupplied Jacobian function, the Jacobian matrix \(J(t,y)\) is zeroed out, so only nonzero elements need to be loaded into Jac.With the default Newton nonlinear solver, each call to the user’s
ARKLsJacFn()
function is preceded by a call to the implicitARKRhsFn()
user function with the same \((t,y)\) arguments. Thus, the Jacobian function can use any auxiliary data that is computed and saved during the evaluation of \(f^I(t,y)\). In the case of a usersupplied or external nonlinear solver, this is also true if the nonlinear system function is evaluated prior to calling the linear solver setup function (see §11.1.4 for more information).If the user’s
ARKLsJacFn
function uses difference quotient approximations, then it may need to access quantities not in the argument list, including the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to theark_mem
structure to theiruser_data
, and then use theARKSodeGet*
functions listed in §2.4.3.10. The unit roundoff can be accessed asSUN_UNIT_ROUNDOFF
, which is defined in the header filesundials_types.h
.dense \(J(t,y)\): A usersupplied dense Jacobian function must load the N by N dense matrix Jac with an approximation to the Jacobian matrix \(J(t,y)\) at the point \((t,y)\). Utility routines and accessor macros for the SUNMATRIX_DENSE module are documented in §9.9.
banded \(J(t,y)\): A usersupplied banded Jacobian function must load the band matrix Jac with the elements of the Jacobian \(J(t,y)\) at the point \((t,y)\). Utility routines and accessor macros for the SUNMATRIX_BAND module are documented in §9.12.
sparse \(J(t,y)\): A usersupplied sparse Jacobian function must load the compressedsparsecolumn (CSC) or compressedsparserow (CSR) matrix Jac with an approximation to the Jacobian matrix \(J(t,y)\) at the point \((t,y)\). Storage for Jac already exists on entry to this function, although the user should ensure that sufficient space is allocated in Jac to hold the nonzero values to be set; if the existing space is insufficient the user may reallocate the data and index arrays as needed. Utility routines and accessor macros for the SUNMATRIX_SPARSE type are documented in §9.14.

typedef int (*ARKLsLinSysFn)(sunrealtype t, N_Vector y, N_Vector fy, SUNMatrix A, SUNMatrix M, sunbooleantype jok, sunbooleantype *jcur, sunrealtype gamma, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
This function computes the linear system matrix \(\mathcal{A}(t,y) = M(t)  \gamma J(t,y)\) (or an approximation to it).
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector, namely the predicted value of \(y(t)\).
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param A:
the output linear system matrix.
 Param M:
the current mass matrix (this input is
NULL
if \(M = I\)). Param jok:
is an input flag indicating whether the Jacobianrelated data needs to be updated. The jok argument provides for the reuse of Jacobian data. When jok =
SUNFALSE
, the Jacobianrelated data should be recomputed from scratch. When jok =SUNTRUE
the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of gamma). A call with jok =SUNTRUE
can only occur after a call with jok =SUNFALSE
. Param jcur:
is a pointer to a flag which should be set to
SUNTRUE
if Jacobian data was recomputed, or set toSUNFALSE
if Jacobian data was not recomputed, but saved data was still reused. Param gamma:
the scalar \(\gamma\) appearing in the Newton system matrix \(\mathcal{A}=M(t)\gamma J(t,y)\).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Param tmp*:
pointers to memory allocated to variables of type
N_Vector
which can be used by an ARKLsLinSysFn as temporary storage or work space. Return:
An ARKLsLinSysFn function should return 0 if successful, a positive value if a recoverable error occurred (in which case ARKODE will attempt to correct, while ARKLS sets last_flag to ARKLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted,
ARKodeEvolve()
returns ARK_LSETUP_FAIL and ARKLS sets last_flag to ARKLS_JACFUNC_UNRECVR).
2.4.4.9. Jacobianvector product
When using a matrixfree linear solver module for the implicit
stage solves (i.e., a NULLvalued SUNMATRIX argument was supplied to
ARKodeSetLinearSolver()
, the user may provide a function
of type ARKLsJacTimesVecFn
in the following form, to compute
matrixvector products \(Jv\). If such a function is not supplied,
the default is a difference quotient approximation to these products.

typedef int (*ARKLsJacTimesVecFn)(N_Vector v, N_Vector Jv, sunrealtype t, N_Vector y, N_Vector fy, void *user_data, N_Vector tmp)
This function computes the product \(Jv\) where \(J(t,y) \approx \dfrac{\partial f^I}{\partial y}(t,y)\) (or an approximation to it).
 Param v:
the vector to multiply.
 Param Jv:
the output vector computed.
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Param tmp:
pointer to memory allocated to a variable of type
N_Vector
which can be used as temporary storage or work space. Return:
The value to be returned by the Jacobianvector product function should be 0 if successful. Any other return value will result in an unrecoverable error of the generic Krylov solver, in which case the integration is halted.
Note
If the user’s
ARKLsJacTimesVecFn
function uses difference quotient approximations, it may need to access quantities not in the argument list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to theark_mem
structure to theiruser_data
, and then use theARKodeGet*
functions listed in §2.4.3.10. The unit roundoff can be accessed asSUN_UNIT_ROUNDOFF
, which is defined in the header filesundials_types.h
.
2.4.4.10. Jacobianvector product setup
If the user’s Jacobiantimesvector routine requires that any Jacobianrelated data
be preprocessed or evaluated, then this needs to be done in a
usersupplied function of type ARKLsJacTimesSetupFn
,
defined as follows:

typedef int (*ARKLsJacTimesSetupFn)(sunrealtype t, N_Vector y, N_Vector fy, void *user_data)
This function preprocesses and/or evaluates any Jacobianrelated data needed by the Jacobiantimesvector routine.
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
The value to be returned by the Jacobianvector setup function should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).
Note
Each call to the Jacobianvector setup function is preceded by a call to the implicit
ARKRhsFn
user function with the same \((t,y)\) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the implicit ODE righthand side.If the user’s
ARKLsJacTimesSetupFn
function uses difference quotient approximations, it may need to access quantities not in the argument list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to theark_mem
structure to theiruser_data
, and then use theARKodeGet*
functions listed in §2.4.3.10. The unit roundoff can be accessed asSUN_UNIT_ROUNDOFF
, which is defined in the header filesundials_types.h
.
2.4.4.11. Preconditioner solve
If a usersupplied preconditioner is to be used with a SUNLinSol
solver module, then the user must provide a function of type
ARKLsPrecSolveFn
to solve the linear system \(Pz=r\),
where \(P\) corresponds to either a left or right
preconditioning matrix. Here \(P\) should approximate (at least
crudely) the Newton matrix \(\mathcal{A}(t,y)=M(t)\gamma J(t,y)\),
where \(M(t)\) is the mass matrix and \(J(t,y) = \dfrac{\partial f^I}{\partial
y}(t,y)\) If preconditioning is done on both sides, the product of the two
preconditioner matrices should approximate \(\mathcal{A}\).

typedef int (*ARKLsPrecSolveFn)(sunrealtype t, N_Vector y, N_Vector fy, N_Vector r, N_Vector z, sunrealtype gamma, sunrealtype delta, int lr, void *user_data)
This function solves the preconditioner system \(Pz=r\).
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param r:
the righthand side vector of the linear system.
 Param z:
the computed output solution vector.
 Param gamma:
the scalar \(\gamma\) appearing in the Newton matrix given by \(\mathcal{A}=M(t)\gamma J(t,y)\).
 Param delta:
an input tolerance to be used if an iterative method is employed in the solution. In that case, the residual vector \(Res = rPz\) of the system should be made to be less than delta in the weighted \(l_2\) norm, i.e. \(\left(\displaystyle \sum_{i=1}^n \left(Res_i * ewt_i\right)^2 \right)^{1/2} < \delta\), where \(\delta =\) delta. To obtain the
N_Vector
ewt, callARKodeGetErrWeights()
. Param lr:
an input flag indicating whether the preconditioner solve is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
The value to be returned by the preconditioner solve function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).
2.4.4.12. Preconditioner setup
If the user’s preconditioner routine above requires that any data be
preprocessed or evaluated, then these actions need to occur within a
usersupplied function of type ARKLsPrecSetupFn
.

typedef int (*ARKLsPrecSetupFn)(sunrealtype t, N_Vector y, N_Vector fy, sunbooleantype jok, sunbooleantype *jcurPtr, sunrealtype gamma, void *user_data)
This function preprocesses and/or evaluates Jacobianrelated data needed by the preconditioner.
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param fy:
the current value of the vector \(f^I(t,y)\).
 Param jok:
is an input flag indicating whether the Jacobianrelated data needs to be updated. The jok argument provides for the reuse of Jacobian data in the preconditioner solve function. When jok =
SUNFALSE
, the Jacobianrelated data should be recomputed from scratch. When jok =SUNTRUE
the Jacobian data, if saved from the previous call to this function, can be reused (with the current value of gamma). A call with jok =SUNTRUE
can only occur after a call with jok =SUNFALSE
. Param jcurPtr:
is a pointer to a flag which should be set to
SUNTRUE
if Jacobian data was recomputed, or set toSUNFALSE
if Jacobian data was not recomputed, but saved data was still reused. Param gamma:
the scalar \(\gamma\) appearing in the Newton matrix given by \(\mathcal{A}=M(t)\gamma J(t,y)\).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
The value to be returned by the preconditioner setup function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).
Note
The operations performed by this function might include forming a crude approximate Jacobian, and performing an LU factorization of the resulting approximation to \(\mathcal{A} = M(t)  \gamma J(t,y)\).
With the default nonlinear solver (the native SUNDIALS Newton method), each call to the preconditioner setup function is preceded by a call to the implicit
ARKRhsFn
user function with the same \((t,y)\) arguments. Thus, the preconditioner setup function can use any auxiliary data that is computed and saved during the evaluation of the implicit ODE righthand side. In the case of a usersupplied or external nonlinear solver, this is also true if the nonlinear system function is evaluated prior to calling the linear solver setup function (see §11.1.4 for more information).This function is not called in advance of every call to the preconditioner solve function, but rather is called only as often as needed to achieve convergence in the Newton iteration.
If the user’s
ARKLsPrecSetupFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current step size, the error weights, etc. To obtain these, the user will need to add a pointer to theark_mem
structure to theiruser_data
, and then use theARKodeGet*
functions listed in §2.4.3.10. The unit roundoff can be accessed asSUN_UNIT_ROUNDOFF
, which is defined in the header filesundials_types.h
.
2.4.4.13. Mass matrix construction
For problems involving a nonidentity mass matrix, if a matrixbased
massmatrix linear solver is used (i.e., a nonNULL SUNMATRIX was
supplied to ARKodeSetMassLinearSolver()
, the
user must provide a function of type ARKLsMassFn
to provide
the mass matrix approximation.

typedef int (*ARKLsMassFn)(sunrealtype t, SUNMatrix M, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
This function computes the mass matrix \(M(t)\) (or an approximation to it).
 Param t:
the current value of the independent variable.
 Param M:
the output mass matrix.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Param tmp1*:
pointers to memory allocated to variables of type
N_Vector
which can be used by an ARKLsMassFn as temporary storage or work space. Return:
An ARKLsMassFn function should return 0 if successful, or a negative value if it failed unrecoverably (in which case the integration is halted,
ARKodeEvolve()
returns ARK_MASSSETUP_FAIL and ARKLS sets last_flag to ARKLS_MASSFUNC_UNRECVR).
Note
Information regarding the structure of the specific
SUNMatrix
structure (e.g.~number of rows, upper/lower bandwidth, sparsity type) may be obtained through using the implementationspecificSUNMatrix
interface functions (see §9 for details).Prior to calling the usersupplied mass matrix function, the mass matrix \(M(t)\) is zeroed out, so only nonzero elements need to be loaded into M.
dense \(M(t)\): A usersupplied dense mass matrix function must load the N by N dense matrix M with an approximation to the mass matrix \(M(t)\). Utility routines and accessor macros for the SUNMATRIX_DENSE module are documented in §9.9.
banded \(M(t)\): A usersupplied banded mass matrix function must load the band matrix M with the elements of the mass matrix \(M(t)\). Utility routines and accessor macros for the SUNMATRIX_BAND module are documented in §9.12.
sparse \(M(t)\): A usersupplied sparse mass matrix function must load the compressedsparsecolumn (CSR) or compressedsparserow (CSR) matrix M with an approximation to the mass matrix \(M(t)\). Storage for M already exists on entry to this function, although the user should ensure that sufficient space is allocated in M to hold the nonzero values to be set; if the existing space is insufficient the user may reallocate the data and row index arrays as needed. Utility routines and accessor macros for the SUNMATRIX_SPARSE type are documented in §9.14.
2.4.4.14. Mass matrixvector product
For problems involving a nonidentity mass matrix, if a matrixfree
linear solver is to be used for massmatrix linear systems (i.e., a
NULLvalued SUNMATRIX argument was supplied to
ARKodeSetMassLinearSolver()
in
§2.4.2), the user must provide a
function of type ARKLsMassTimesVecFn
in the following form, to
compute matrixvector products \(M(t)\, v\).

typedef int (*ARKLsMassTimesVecFn)(N_Vector v, N_Vector Mv, sunrealtype t, void *mtimes_data)
This function computes the product \(M(t)\, v\) (or an approximation to it).
 Param v:
the vector to multiply.
 Param Mv:
the output vector computed.
 Param t:
the current value of the independent variable.
 Param mtimes_data:
a pointer to user data, the same as the mtimes_data parameter that was passed to
ARKodeSetMassTimes()
. Return:
The value to be returned by the massmatrixvector product function should be 0 if successful. Any other return value will result in an unrecoverable error of the generic Krylov solver, in which case the integration is halted.
2.4.4.15. Mass matrixvector product setup
For problems involving a nonidentity mass matrix and a matrixfree linear
solver, if the user’s massmatrixtimesvector routine requires that any
mass matrixrelated data be preprocessed or evaluated, then this needs to
be done in a usersupplied function of type
ARKLsMassTimesSetupFn
, defined as follows:

typedef int (*ARKLsMassTimesSetupFn)(sunrealtype t, void *mtimes_data)
This function preprocesses and/or evaluates any massmatrixrelated data needed by the massmatrixtimesvector routine.
 Param t:
the current value of the independent variable.
 Param mtimes_data:
a pointer to user data, the same as the mtimes_data parameter that was passed to
ARKodeSetMassTimes()
. Return:
The value to be returned by the massmatrixvector setup function should be 0 if successful. Any other return value will result in an unrecoverable error of the ARKLS mass matrix solver interface, in which case the integration is halted.
2.4.4.16. Mass matrix preconditioner solve
For problems involving a nonidentity mass matrix and an iterative linear
solver, if a usersupplied preconditioner is to be used with a SUNLINEAR
solver module for mass matrix linear systems, then the user must
provide a function of type ARKLsMassPrecSolveFn
to solve the
linear system \(Pz=r\), where \(P\) may be either a left or right
preconditioning matrix. Here \(P\) should approximate (at least
crudely) the mass matrix \(M(t)\). If preconditioning is done on
both sides, the product of the two preconditioner matrices should
approximate \(M(t)\).

typedef int (*ARKLsMassPrecSolveFn)(sunrealtype t, N_Vector r, N_Vector z, sunrealtype delta, int lr, void *user_data)
This function solves the preconditioner system \(Pz=r\).
 Param t:
the current value of the independent variable.
 Param r:
the righthand side vector of the linear system.
 Param z:
the computed output solution vector.
 Param delta:
an input tolerance to be used if an iterative method is employed in the solution. In that case, the residual vector \(Res = rPz\) of the system should be made to be less than delta in the weighted \(l_2\) norm, i.e. \(\left(\displaystyle \sum_{i=1}^n \left(Res_i * ewt_i\right)^2 \right)^{1/2} < \delta\), where \(\delta =\) delta. To obtain the
N_Vector
ewt, callARKodeGetErrWeights()
. Param lr:
an input flag indicating whether the preconditioner solve is to use the left preconditioner (lr = 1) or the right preconditioner (lr = 2).
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
The value to be returned by the preconditioner solve function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).
2.4.4.17. Mass matrix preconditioner setup
For problems involving a nonidentity mass matrix and an iterative linear
solver, if the user’s mass matrix preconditioner above requires that any
problem data be preprocessed or evaluated, then these actions need to
occur within a usersupplied function of type
ARKLsMassPrecSetupFn
.

typedef int (*ARKLsMassPrecSetupFn)(sunrealtype t, void *user_data)
This function preprocesses and/or evaluates massmatrixrelated data needed by the preconditioner.
 Param t:
the current value of the independent variable.
 Param user_data:
a pointer to user data, the same as the user_data parameter that was passed to
ARKodeSetUserData()
. Return:
The value to be returned by the mass matrix preconditioner setup function is a flag indicating whether it was successful. This value should be 0 if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).
Note
The operations performed by this function might include forming a mass matrix and performing an incomplete factorization of the result. Although such operations would typically be performed only once at the beginning of a simulation, these may be required if the mass matrix can change as a function of time.
If both this function and a
ARKLsMassTimesSetupFn
are supplied, all calls to this function will be preceded by a call to theARKLsMassTimesSetupFn
, so any setup performed there may be reused.
2.4.4.18. Vector resize function
For simulations involving changes to the number of equations and
unknowns in the ODE system (e.g. when using spatial adaptivity in a
PDE simulation), the ARKODE integrator may be “resized” between
integration steps, through calls to the ARKodeResize()
function. Typically, when performing adaptive simulations the solution
is stored in a customized usersupplied data structure, to enable
adaptivity without repeated allocation/deallocation of memory. In
these scenarios, it is recommended that the user supply a customized
vector kernel to interface between SUNDIALS and their problemspecific
data structure. If this vector kernel includes a function of type
ARKVecResizeFn
to resize a given vector implementation, then
this function may be supplied to ARKodeResize()
so that all
internal ARKODE vectors may be resized, instead of deleting and
recreating them at each call. This resize function should have the
following form:

typedef int (*ARKVecResizeFn)(N_Vector y, N_Vector ytemplate, void *user_data)
This function resizes the vector y to match the dimensions of the supplied vector, ytemplate.
 Param y:
the vector to resize.
 Param ytemplate:
a vector of the desired size.
 Param user_data:
a pointer to user data, the same as the resize_data parameter that was passed to
ARKodeResize()
. Return:
An ARKVecResizeFn function should return 0 if it successfully resizes the vector y, and a nonzero value otherwise.
Note
If this function is not supplied, then ARKODE will instead destroy the vector y and clone a new vector y off of ytemplate.
2.4.4.19. Pre inner integrator communication function (MRIStep only)
The user may supply a function of type MRIStepPreInnerFn
that will be
called before each inner integration to perform any communication or
memory transfers of forcing data supplied by the outer integrator to the inner
integrator for the inner integration.

typedef int (*MRIStepPreInnerFn)(sunrealtype t, N_Vector *f, int num_vecs, void *user_data)
 Param t:
the current value of the independent variable.
 Param f:
an
N_Vector
array of outer forcing vectors. Param num_vecs:
the number of vectors in the
N_Vector
array. Param user_data:
the user_data pointer that was passed to
MRIStepSetUserData()
. Return:
An MRIStepPreInnerFn function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. As the MRIStep module only supports fixed step sizes at this time any nonzero return value will halt the integration.
Note
In a heterogeneous computing environment if any data copies between the host and device vector data are necessary, this is where that should occur.
2.4.4.20. Post inner integrator communication function (MRIStep only)
The user may supply a function of type MRIStepPostInnerFn
that will be
called after each inner integration to perform any communication or
memory transfers of state data supplied by the inner integrator to the
outer integrator for the outer integration.

typedef int (*MRIStepPostInnerFn)(sunrealtype t, N_Vector y, void *user_data)
 Param t:
the current value of the independent variable.
 Param y:
the current value of the dependent variable vector.
 Param user_data:
the
user_data
pointer that was passed toMRIStepSetUserData()
. Return:
An
MRIStepPostInnerFn()
function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. As the MRIStep module only supports fixed step sizes at this time any nonzero return value will halt the integration.
Note
In a heterogeneous computing environment if any data copies between the host and device vector data are necessary, this is where that should occur.
2.4.4.21. Relaxation function

typedef int (*ARKRelaxFn)(N_Vector y, sunrealtype *r, void *user_data)
When applying relaxation, an
ARKRelaxFn()
function is required to compute the conservative or dissipative function \(\xi(y)\). Param y:
the current value of the dependent variable vector.
 Param r:
the value of \(\xi(y)\).
 Param user_data:
the
user_data
pointer that was passed toARKodeSetUserData()
. Return:
An
ARKRelaxFn()
function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. If a recoverable error occurs, the step size will be reduced and the step repeated.
2.4.4.22. Relaxation Jacobian function

typedef int (*ARKRelaxJacFn)(N_Vector y, N_Vector J, void *user_data);
When applying relaxation, an
ARKRelaxJacFn()
function is required to compute the Jacobian \(\xi'(y)\) of theARKRelaxFn()
\(\xi(y)\). Param y:
the current value of the dependent variable vector.
 Param J:
the Jacobian vector \(\xi'(y)\).
 Param user_data:
the
user_data
pointer that was passed toARKodeSetUserData()
. Return:
An
ARKRelaxJacFn()
function should return 0 if successful, a positive value if a recoverable error occurred, or a negative value if an unrecoverable error occurred. If a recoverable error occurs, the step size will be reduced and the step repeated.