6.4. Using IDA for IVP Solution
This chapter is concerned with the use of IDA for the integration of DAEs.
The following sections treat the header files and the layout of the user’s main program, and provide descriptions of the IDA usercallable functions and usersupplied functions. The sample programs described in the companion document [69] may also be helpful. Those codes may be used as templates (with the removal of some lines used in testing) and are included in the IDA package.
IDA uses various constants for both input and output. These are defined as needed in this chapter, but for convenience are also listed separately in §6.5.
The user should be aware that not all SUNLinearSolver
and SUNMatrix
objects are compatible with all N_Vector
implementations. Details on
compatibility are given in the documentation for each SUNMatrix
(Chapter
§10) and SUNLinearSolver
(Chapter §11)
implementation. For example, NVECTOR_PARALLEL
is not compatible with the
dense, banded, or sparse SUNMatrix
types, or with the corresponding dense,
banded, or sparse SUNLinearSolver
objects. Please check Chapters
§10 and §11 to verify compatibility between
these objects. In addition to that documentation, we note that the IDABBDPRE
preconditioner can only be used with NVECTOR_PARALLEL
. It is not recommended
to use a threaded vector object with SuperLU_MT unless it is the
NVECTOR_OPENMP
module, and SuperLU_MT is also compiled with OpenMP.
6.4.1. Access to library and header files
At this point, it is assumed that the installation of IDA, following the procedure described in §14, has been completed successfully.
Regardless of where the user’s application program resides, its associated compilation and load commands must make reference to the appropriate locations for the library and header files required by IDA. The relevant library files are
<libdir>/libsundials_ida.<soa>
<libdir>/libsundials_nvec*.<soa>
<libdir>/libsundials_sunmat*.<soa>
<libdir>/libsundials_sunlinsol*.<soa>
<libdir>/libsundials_sunnonlinsol*.<soa>
where the file extension .so
is typically for shared libraries and .a
for static libraries. The relevant header files are located in the
subdirectories
<incdir>/ida
<incdir>/sundials
<incdir>/nvector
<incdir>/sunmatrix
<incdir>/sunlinsol
<incdir>/sunnonlinsol
The directories libdir
and incdir
are the install library and include
directories, respectively. For a default installation, these are
<instdir>/lib
or <instdir>/lib64
and <instdir>/include
,
respectively, where instdir
is the directory where SUNDIALS was installed
(see §14).
6.4.2. Data Types
The header file sundials_types.h
contains the definition of the types:
realtype
– the floatingpoint type used by the SUNDIALS packagessunindextype
– the integer type used for vector and matrix indicesbooleantype
– the type used for logic operations within SUNDIALSSUNOutputFormat
– an enumerated type for SUNDIALS output formats
6.4.2.1. Floating point types

type realtype
The type
realtype
can befloat
,double
, orlong double
, with the default beingdouble
. The user can change the precision of the arithmetic used in the SUNDIALS solvers at the configuration stage (seeSUNDIALS_PRECISION
).
Additionally, based on the current precision, sundials_types.h
defines
BIG_REAL
to be the largest value representable as a realtype
,
SMALL_REAL
to be the smallest value representable as a realtype
, and
UNIT_ROUNDOFF
to be the difference between \(1.0\) and the minimum
realtype
greater than \(1.0\).
Within SUNDIALS, real constants are set by way of a macro called RCONST
. It
is this macro that needs the ability to branch on the definition of
realtype
. In ANSI C, a floatingpoint constant with no suffix is stored as a
double
. Placing the suffix “F
” at the end of a floating point constant
makes it a float
, whereas using the suffix “L
” makes it a long
double
. For example,
#define A 1.0
#define B 1.0F
#define C 1.0L
defines A
to be a double
constant equal to \(1.0\), B
to be a
float
constant equal to \(1.0\), and C
to be a long double
constant equal to \(1.0\). The macro call RCONST(1.0)
automatically
expands to 1.0
if realtype
is double
, to 1.0F
if realtype
is
float
, or to 1.0L
if realtype
is long double
. SUNDIALS uses the
RCONST
macro internally to declare all of its floatingpoint constants.
Additionally, SUNDIALS defines several macros for common mathematical functions
e.g., fabs
, sqrt
, exp
, etc. in sundials_math.h
. The macros are
prefixed with SUNR
and expand to the appropriate C
function based on the
realtype
. For example, the macro SUNRabs
expands to the C
function
fabs
when realtype
is double
, fabsf
when realtype
is
float
, and fabsl
when realtype
is long double
.
A user program which uses the type realtype
, the RCONST
macro, and the
SUNR
mathematical function macros is precisionindependent except for any
calls to precisionspecific library functions. Our example programs use
realtype
, RCONST
, and the SUNR
macros. Users can, however, use the
type double
, float
, or long double
in their code (assuming that this
usage is consistent with the typedef for realtype
) and call the appropriate
math library functions directly. Thus, a previously existing piece of C or C++
code can use SUNDIALS without modifying the code to use realtype
,
RCONST
, or the SUNR
macros so long as the SUNDIALS libraries are built
to use the corresponding precision (see §14.1.2).
6.4.2.2. Integer types used for indexing

type sunindextype
The type
sunindextype
is used for indexing array entries in SUNDIALS modules as well as for storing the total problem size (e.g., vector lengths and matrix sizes). During configurationsunindextype
may be selected to be either a 32 or 64bit signed integer with the default being 64bit (seeSUNDIALS_INDEX_SIZE
).
When using a 32bit integer the total problem size is limited to
\(2^{31}1\) and with 64bit integers the limit is \(2^{63}1\). For
users with problem sizes that exceed the 64bit limit an advanced configuration
option is available to specify the type used for sunindextype
(see SUNDIALS_INDEX_TYPE
).
A user program which uses sunindextype
to handle indices will work with both
index storage types except for any calls to index storagespecific external
libraries. Our C
and C++
example programs use sunindextype
. Users
can, however, use any compatible type (e.g., int
, long int
,
int32_t
, int64_t
, or long long int
) in their code, assuming that
this usage is consistent with the typedef for sunindextype
on their
architecture. Thus, a previously existing piece of C or C++ code can use
SUNDIALS without modifying the code to use sunindextype
, so long as the
SUNDIALS libraries use the appropriate index storage type (for details see
§14.1.2).
6.4.2.3. Boolean type

type booleantype
As ANSI C89 (ISO C90) does not have a builtin boolean data type, SUNDIALS defines the type
booleantype
as anint
.
The advantage of using the name booleantype (instead of int) is an increase in
code readability. It also allows the programmer to make a distinction between
int and boolean data. Variables of type booleantype
are intended to have
only the two values SUNFALSE
(0
) and SUNTRUE
(1
).
6.4.2.4. Output formatting type

enum SUNOutputFormat
The enumerated type
SUNOutputFormat
defines the enumeration constants for SUNDIALS output formats

enumerator SUN_OUTPUTFORMAT_TABLE
The output will be a table of values

enumerator SUN_OUTPUTFORMAT_CSV
The output will be a commaseparated list of key and value pairs e.g.,
key1,value1,key2,value2,...
Note
The file
scripts/sundials_csv.py
provides python utility functions to read and output the data from a SUNDIALS CSV output file using the key and value pair format.
6.4.3. Header files
The calling program must include several header files so that various macros and data types can be used. The header file that is always required is:
ida/ida.h
the main header file for IDA, which defines the types and various constants, and includes function prototypes. This includes the header file for IDALS,ida/ida_ls.h
.
Note that ida.h
includes sundials_types.h
, which defines the types,
realtype
, sunindextype
, and booleantype
and the constants
SUNFALSE
and SUNTRUE
.
The calling program must also include an N_Vector
implementation
header file, of the form nvector/nvector_*.h
(see Chapter §9
for more information). This file in turn includes the header file
sundials_nvector.h
which defines the abstract vector data type.
If using a nondefault nonlinear solver object, or when interacting with a
SUNNonlinearSolver
object directly, the calling program must also include a
SUNNonlinearSolver
implementation header file, of the form
sunnonlinsol/sunnonlinsol_*.h
where *
is the name of the nonlinear
solver (see Chapter §12 for more information). This file in
turn includes the header file sundials_nonlinearsolver.h
which defines the
abstract nonlinear linear solver data type.
If using a nonlinear solver that requires the solution of a linear system of the
form (6.3) (e.g., the default Newton iteration), the calling program
must also include a SUNLinearSolver
implementation header file, of the from
sunlinsol/sunlinsol_*.h
where *
is the name of the linear solver
(see Chapter §11 for more information). This file in
turn includes the header file sundials_linearsolver.h
which defines the
abstract linear solver data type.
If the linear solver is matrixbased, the linear solver header will also include
a header file of the from sunmatrix/sunmatrix_*.h
where *
is the name of
the matrix implementation compatible with the linear solver. The matrix header
file provides access to the relevant matrix functions/macros and in turn
includes the header file sundials_matrix.h
which defines the abstract matrix
data type.
Other headers may be needed, according to the choice of preconditioner, etc. For
example, in the example idaFoodWeb_kry_p
(see [69]),
preconditioning is done with a blockdiagonal matrix. For this, even though the
SUNLINSOL_SPGMR
linear solver is used, the header
sundials/sundials_dense.h
is included for access to the underlying generic
dense matrix arithmetic routines.
6.4.4. A skeleton of the user’s main program
The following is a skeleton of the user’s main program (or calling program) for
the integration of a DAE IVP. Most of the steps are independent of the
N_Vector
, SUNMatrix
, SUNLinearSolver
, and
SUNNonlinearSolver
implementations used. For the steps that are not,
refer to Chapters §9, §10, §11,
and §12 for the specific name of the function to be called or
macro to be referenced.
Initialize parallel or multithreaded environment (if appropriate)
For example, call
MPI_Init
to initialize MPI if used.Create the SUNDIALS context object
Call
SUNContext_Create()
to allocate theSUNContext
object.Create the vector of initial values
Construct an
N_Vector
of initial values using the appropriate functions defined by the particularN_Vector
implementation (see §9 for details).For native SUNDIALS vector implementations, use a call of the form
y0 = N_VMake_***(..., ydata)
if the array containing the initial values of \(y\) already exists. Otherwise, create a new vector by making a call of the formN_VNew_***(...)
, and then set its elements by accessing the underlying data with a call of the formydata = N_VGetArrayPointer(y0)
. Here,***
is the name of the vector implementation.For hypre, PETSc, and Trilinos vector wrappers, first create and initialize the underlying vector, and then create an
N_Vector
wrapper with a call of the formy0 = N_VMake_***(yvec)
, whereyvec
is a hypre, PETSc, or Trilinos vector. Note that calls likeN_VNew_***(...)
andN_VGetArrayPointer(...)
are not available for these vector wrappers.Set the vector
yp0
of initial conditions for \(\dot{y}\) similarly.Create matrix object (if appropriate)
If a linear solver is required (e.g., when using the default Newton solver) and the linear solver will be a matrixbased linear solver, then a template Jacobian matrix must be created by calling the appropriate constructor defined by the particular
SUNMatrix
implementation.For the native SUNDIALS
SUNMatrix
implementations, the matrix object may be created using a call of the formSUN***Matrix(...)
where***
is the name of the matrix (see §10 for details).Create linear solver object (if appropriate)
If a linear solver is required (e.g., when using the default Newton solver), then the desired linear solver object must be created by calling the appropriate constructor defined by the particular
SUNLinearSolver
implementation.For any of the native SUNDIALS
SUNLinearSolver
implementations, the linear solver object may be created using a call of the formSUNLinearSolver LS = SUNLinSol_***(...);
where***
is the name of the linear solver (see §11 for details).Create nonlinear solver object (if appropriate)
If using a nondefault nonlinear solver, then the desired nonlinear solver object must be created by calling the appropriate constructor defined by the particular
SUNNonlinearSolver
implementation.For any of the native SUNDIALS
SUNNonLinearSolver
implementations, the nonlinear solver object may be created using a call of the formSUNNonlinearSolver NLS = SUNNonlinSol_***(...);
where***
is the name of the nonlinear solver (see §12 for details).Create IDA object
Call
IDACreate()
to create the IDA solver object.Initialize IDA solver
Call
IDAInit()
to provide the initial condition vectors created above, set the DAE residual function, and initialize IDA.Specify integration tolerances
Call one of the following functions to set the integration tolerances:
IDASStolerances()
to specify scalar relative and absolute tolerances.IDASVtolerances()
to specify a scalar relative tolerance and a vector of absolute tolerances.IDAWFtolerances()
to specify a function which sets directly the weights used in evaluating WRMS vector norms.
See §6.4.5.3 for general advice on selecting tolerances and §6.4.5.4 for advice on controlling unphysical values.
Attach the linear solver (if appropriate)
If a linear solver was created above, initialize the IDALS linear solver interface by attaching the linear solver object (and matrix object, if applicable) with
IDASetLinearSolver()
.Set linear solver optional inputs (if appropriate)
See Table 6.2 for IDALS optional inputs and Chapter §11 for linear solver specific optional inputs.
Attach nonlinear solver module (if appropriate)
If a nonlinear solver was created above, initialize the IDANLS nonlinear solver interface by attaching the nonlinear solver object with
IDASetNonlinearSolver()
.Set nonlinear solver optional inputs (if appropriate)
See Table 6.3 for IDANLS optional inputs and Chapter §12 for nonlinear solver specific optional inputs. Note, solver specific optional inputs must be called after
IDASetNonlinearSolver()
, otherwise the optional inputs will be overridden by IDA defaults.Specify rootfinding problem (optional)
Call
IDARootInit()
to initialize a rootfinding problem to be solved during the integration of the ODE system. See Table 6.6 for relevant optional input calls.Set optional inputs
Call
IDASet***
functions to change any optional inputs that control the behavior of IDA from their default values. See §6.4.5.10 for details.Correct initial values (optional)
Call
IDACalcIC()
to correct the initial valuesy0
andyp0
passed toIDAInit()
. See Table 6.4 for relevant optional input calls.Advance solution in time
For each point at which output is desired, call
ier = IDASolve(ida_mem, tout, &tret, yret, ypret, itask)
. Hereitask
specifies the return mode. The vectoryret
(which can be the same as the vectory0
above) will contain \(y(t)\), while the vectorypret
(which can be the same as the vectoryp0
above) will contain \(\dot{y}(t)\).See
IDASolve()
for details.Get optional outputs
Call
IDAGet***
functions to obtain optional output. See §6.4.5.12 for details.Destroy objects
Upon completion of the integration call the following functions, as necessary, to destroy any objects created above:
Call
N_VDestroy()
to free vector objects.Call
SUNMatDestroy()
to free matrix objects.Call
SUNLinSolFree()
to free linear solvers objects.Call
SUNNonlinSolFree()
to free nonlinear solvers objects.Call
IDAFree()
to free the memory allocated by IDA.Call
SUNContext_Free()
to free the SUNDIALS context.
Finalize MPI, if used
Call
MPI_Finalize
to terminate MPI.
6.4.5. Usercallable functions
This section describes the IDA functions that are called by the user to setup and then solve an IVP. Some of these are required. However, starting with §6.4.5.10, the functions listed involve optional inputs/outputs or restarting, and those paragraphs may be skipped for a casual use of IDA. In any case, refer to §6.4.4 for the correct order of these calls.
On an error, each usercallable function returns a negative value and sends an
error message to the error handler routine, which prints the message on
stderr
by default. However, the user can set a file as error output or can
provide his own error handler function (see
§6.4.5.10.1).
6.4.5.1. IDA initialization and deallocation functions

void *IDACreate(SUNContext sunctx)
The function
IDACreate
instantiates an IDA solver object. Arguments:
sunctx
– theSUNContext
object (see §2.1)
 Return value:
void*
pointer the IDA solver object.

int IDAInit(void *ida_mem, IDAResFn res, realtype t0, N_Vector y0, N_Vector yp0)
The function
IDAInit
provides required problem and solution specifications, allocates internal memory, and initializes IDA. Arguments:
ida_mem
– pointer to the IDA solver object.res
– is the function which computes the residual function \(F(t, y, \dot{y})\) for the DAE. For full details seeIDAResFn
.t0
– is the initial value of \(t\).y0
– is the initial value of \(y\).yp0
– is the initial value of \(\dot{y}\).
 Return value:
IDA_SUCCESS
– The call was successful.IDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_MEM_FAIL
– A memory allocation request has failed.IDA_ILL_INPUT
– An input argument toIDAInit()
has an illegal value.
 Notes:
If an error occurred,
IDAInit()
also sends an error message to the error handler function.

void IDAFree(void **ida_mem)
The function
IDAFree
frees the pointer allocated by a previous call toIDACreate()
. Arguments:
ida_mem
– pointer to the IDA solver object.
 Return value:
void
6.4.5.2. IDA tolerance specification functions
One of the following three functions must be called to specify the integration
tolerances (or directly specify the weights used in evaluating WRMS vector
norms). Note that this call must be made after the call to IDAInit()
.

int IDASStolerances(void *ida_mem, realtype reltol, realtype abstol)
The function
IDASStolerances
specifies scalar relative and absolute tolerances. Arguments:
ida_mem
– pointer to the IDA solver object.reltol
– is the scalar relative error tolerance.abstol
– is the scalar absolute error tolerance.
 Return value:
IDA_SUCCESS
– The call was successfulIDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_NO_MALLOC
– The allocation functionIDAInit()
has not been called.IDA_ILL_INPUT
– One of the input tolerances was negative.

int IDASVtolerances(void *ida_mem, realtype reltol, N_Vector abstol)
The function
IDASVtolerances
specifies scalar relative tolerance and vector absolute tolerances. Arguments:
ida_mem
– pointer to the IDA solver object.reltol
– is the scalar relative error tolerance.abstol
– is the vector of absolute error tolerances.
 Return value:
IDA_SUCCESS
– The call was successfulIDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_NO_MALLOC
– The allocation functionIDAInit()
has not been called.IDA_ILL_INPUT
– The relative error tolerance was negative or the absolute tolerance vector had a negative component.
 Notes:
This choice of tolerances is important when the absolute error tolerance needs to be different for each component of the state vector \(y\).

int IDAWFtolerances(void *ida_mem, IDAEwtFn efun)
The function
IDAWFtolerances
specifies a usersupplied functionefun
that sets the multiplicative error weights \(W_i\) for use in the weighted RMS norm, which are normally defined by (6.4). Arguments:
ida_mem
– pointer to the IDA solver object.IDACreate()
efun
– is the function which defines theewt
vector. For full details seeIDAEwtFn
.
 Return value:
IDA_SUCCESS
– The call was successfulIDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_NO_MALLOC
– The allocation functionIDAInit()
has not been called.
6.4.5.3. General advice on choice of tolerances
For many users, the appropriate choices for tolerance values in reltol
and
abstol
are a concern. The following pieces of advice are relevant.
The scalar relative tolerance
reltol
is to be set to control relative errors. Soreltol
of \(10^{4}\) means that errors are controlled to .01%. We do not recommend usingreltol
larger than \(10^{3}\). On the other hand,reltol
should not be so small that it is comparable to the unit roundoff of the machine arithmetic (generally around \(10^{15}\)).The absolute tolerances
abstol
(whether scalar or vector) need to be set to control absolute errors when any components of the solution vectory
may be so small that pure relative error control is meaningless. For example, ify[i]
starts at some nonzero value, but in time decays to zero, then pure relative error control ony[i]
makes no sense (and is overly costly) aftery[i]
is below some noise level. Thenabstol
(if a scalar) orabstol[i]
(if a vector) needs to be set to that noise level. If the different components have different noise levels, thenabstol
should be a vector. See the exampleidaRoberts_dns
in the IDA package, and the discussion of it in the IDA Examples document [69]. In that problem, the three components vary betwen 0 and 1, and have different noise levels; hence theabstol
vector. It is impossible to give any general advice onabstol
values, because the appropriate noise levels are completely problemdependent. The user or modeler hopefully has some idea as to what those noise levels are.Finally, it is important to pick all the tolerance values conservatively, because they control the error committed on each individual time step. The final (global) errors are some sort of accumulation of those perstep errors. A good rule of thumb is to reduce the tolerances by a factor of .01 from the actual desired limits on errors. So if you want .01% accuracy (globally), a good choice is to is a
reltol
of \(10^{6}\). But in any case, it is a good idea to do a few experiments with the tolerances to see how the computed solution values vary as tolerances are reduced.
6.4.5.4. Advice on controlling unphysical negative values
In many applications, some components in the true solution are always positive or nonnegative, though at times very small. In the numerical solution, however, small negative (hence unphysical) values can then occur. In most cases, these values are harmless, and simply need to be controlled, not eliminated. The following pieces of advice are relevant.
The way to control the size of unwanted negative computed values is with tighter absolute tolerances. Again this requires some knowledge of the noise level of these components, which may or may not be different for different components. Some experimentation may be needed.
If output plots or tables are being generated, and it is important to avoid having negative numbers appear there (for the sake of avoiding a long explanation of them, if nothing else), then eliminate them, but only in the context of the output medium. Then the internal values carried by the solver are unaffected. Remember that a small negative value in
yret
returned by IDA, with magnitude comparable toabstol
or less, is equivalent to zero as far as the computation is concerned.The user’s residual function
res
should never change a negative value in the solution vectoryy
to a nonnegative value, as a “solution” to this problem. This can cause instability. If theres
routine cannot tolerate a zero or negative value (e.g., because there is a square root or log of it), then the offending value should be changed to zero or a tiny positive number in a temporary variable (not in the inputyy
vector) for the purposes of computing \(F(t,y,\dot{y})\).IDA provides the option of enforcing positivity or nonnegativity on components. Also, such constraints can be enforced by use of the recoverable error return feature in the usersupplied residual function. However, because these options involve some extra overhead cost, they should only be exercised if the use of absolute tolerances to control the computed values is unsuccessful.
6.4.5.5. Linear solver interface functions
As previously explained, if the nonlinear solver requires the solution of linear
systems of the form (6.5), e.g., the default Newton solver, then
the solution of these linear systems is handled with the IDALS linear solver
interface. This interface supports all valid SUNLinearSolver
objects.
Here, a matrixbased SUNLinearSolver
utilizes SUNMatrix
objects to store the Jacobian matrix \(J = \dfrac{\partial{F}}{\partial{y}} + \alpha
\dfrac{\partial{F}}{\partial{\dot{y}}}\) and factorizations used throughout the solution
process. Conversely, matrixfree SUNLinearSolver
object instead use
iterative methods to solve the linear systems of equations, and only require the
action of the Jacobian on a vector, \(Jv\).
With most iterative linear solvers, preconditioning can be done on the left only, on the right only, on both the left and the right, or not at all. The exceptions to this rule are SPFGMR that supports right preconditioning only and PCG that performs symmetric preconditioning. However, in IDA only left preconditioning is supported. For the specification of a preconditioner, see the iterative linear solver sections in §6.4.5.10 and §6.4.6. A preconditioner matrix \(P\) must approximate the Jacobian \(J\), at least crudely.
To attach a generic linear solver to IDA, after the call to IDACreate()
but before any calls to IDASolve()
, the user’s program must create the
appropriate SUNLinearSolver
object and call the function
IDASetLinearSolver()
. To create the SUNLinearSolver
object,
the user may call one of the SUNDIALSpackaged SUNLinearSolver
constructors via a call of the form
SUNLinearSolver LS = SUNLinSol_*(...);
Alternately, a usersupplied SUNLinearSolver
object may be created and
used instead. The use of each of the generic linear solvers involves certain
constants, functions and possibly some macros, that are likely to be needed in
the user code. These are available in the corresponding header file associated
with the specific SUNMatrix
or SUNLinearSolver
object in
question, as described in Chapters §10 and §11.
Once this solver object has been constructed, the user should attach it to IDA
via a call to IDASetLinearSolver()
. The first argument passed to this
function is the IDA memory pointer returned by IDACreate()
; the second
argument is the desired SUNLinearSolver
object to use for solving
systems. The third argument is an optional SUNMatrix
object to
accompany matrixbased SUNLinearSolver
inputs (for matrixfree linear
solvers, the third argument should be NULL
). A call to this function
initializes the IDALS linear solver interface, linking it to the main IDA
integrator, and allows the user to specify additional parameters and routines
pertinent to their choice of linear solver.

int IDASetLinearSolver(void *ida_mem, SUNLinearSolver LS, SUNMatrix J)
The function
IDASetLinearSolver
attaches aSUNLinearSolver
objectLS
and corresponding template JacobianSUNMatrix
objectJ
(if applicable) to IDA, initializing the IDALS linear solver interface. Arguments:
ida_mem
– pointer to the IDA solver object.LS
–SUNLinearSolver
object to use for solving linear systems of the form (6.5).J
–SUNMatrix
object for used as a template for the Jacobian orNULL
if not applicable.
 Return value:
IDALS_SUCCESS
– The IDALS initialization was successful.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_ILL_INPUT
– The IDALS interface is not compatible with theLS
orJ
input objects or is incompatible with theN_Vector
object passed toIDAInit()
.IDALS_SUNLS_FAIL
– A call to theLS
object failed.IDALS_MEM_FAIL
– A memory allocation request failed.
 Notes:
If
LS
is a matrixbased linear solver, then the template Jacobian matrixJ
will be used in the solve process, so if additional storage is required within theSUNMatrix
object (e.g., for factorization of a banded matrix), ensure that the input object is allocated with sufficient size (see the documentation of the particularSUNMatrix
in Chapter §10 for further information).
Warning
The previous routines
IDADlsSetLinearSolver
andIDASpilsSetLinearSolver
are now wrappers for this routine, and may still be used for backwardcompatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
6.4.5.6. Nonlinear solver interface function
By default IDA uses the SUNNonlinearSolver
implementation of Newton’s method
(see §12.7). To attach a different nonlinear solver in
IDA, the user’s program must create a SUNNonlinearSolver
object by calling
the appropriate constructor routine. The user must then attach the
SUNNonlinearSolver
object to IDA by calling IDASetNonlinearSolver()
.
When changing the nonlinear solver in IDA, IDASetNonlinearSolver()
must
be called after IDAInit()
. If any calls to IDASolve()
have been
made, then IDA will need to be reinitialized by calling IDAReInit()
to
ensure that the nonlinear solver is initialized correctly before any subsequent
calls to IDASolve()
.
The first argument passed to IDASetNonlinearSolver()
is the IDA memory
pointer returned by IDACreate()
and the second argument is the
SUNNonlinearSolver
object to use for solving the nonlinear system
(6.3). A call to this function attaches the nonlinear solver to the main
IDA integrator. We note that at present, the SUNNonlinearSolver
object
must be of type SUNNONLINEARSOLVER_ROOTFIND
.

int IDASetNonlinearSolver(void *ida_mem, SUNNonlinearSolver NLS)
The function
IDASetNonLinearSolver
attaches aSUNNonlinearSolver
object (NLS
) to IDA. Arguments:
ida_mem
– pointer to the IDA solver object.NLS
–SUNNonlinearSolver
object to use for solving nonlinear systems.
 Return value:
IDA_SUCCESS
– The nonlinear solver was successfully attached.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– TheSUNNonlinearSolver
object isNULL
, does not implement the required nonlinear solver operations, is not of the correct type, or the residual function, convergence test function, or maximum number of nonlinear iterations could not be set.
6.4.5.7. Initial condition calculation function
IDACalcIC()
calculates corrected initial conditions for the DAE system
for certain indexone problems including a class of systems of semiimplicit
form (see §6.2.2 and [24]). It uses a Newton
iteration combined with a linesearch algorithm. Calling IDACalcIC()
is
optional. It is only necessary when the initial conditions do not satisfy the
given system. Thus if y0
and yp0
are known to satisfy
\(F(t_0, y_0, \dot{y}_0) = 0\), then a call to IDACalcIC()
is
generally not necessary.
A call to the function IDACalcIC()
must be preceded by successful calls
to IDACreate()
and IDAInit()
(or IDAReInit()
), and by a
successful call to the linear system solver specification function. The call to
IDACalcIC()
should precede the call(s) to IDASolve()
for the
given problem.

int IDACalcIC(void *ida_mem, int icopt, realtype tout1)
The function
IDACalcIC
corrects the initial valuesy0
andyp0
at timet0
. Arguments:
ida_mem
– pointer to the IDA solver object.icopt
– is one of the following two options for the initial condition calculation.IDA_YA_YDP_INIT
directsIDACalcIC()
to compute the algebraic components of \(y\) and differential components of \(\dot{y}\), given the differential components of \(y\). This option requires that theN_Vector id
was set throughIDASetId()
, specifying the differential and algebraic components.IDA_Y_INIT
directsIDACalcIC()
to compute all components of \(y\), given \(\dot{y}\). In this case,id
is not required.
tout1
– is the first value of \(t\) at which a solution will be requested (fromIDASolve()
). This value is needed here only to determine the direction of integration and rough scale in the independent variable \(t\).
 Return value:
IDA_SUCCESS
–IDACalcIC()
succeeded.IDA_MEM_NULL
– The argumentida_mem
wasNULL
.IDA_NO_MALLOC
– The allocation functionIDAInit()
has not been called.IDA_ILL_INPUT
– One of the input arguments was illegal.IDA_LSETUP_FAIL
– The linear solver’s setup function failed in an unrecoverable manner.IDA_LINIT_FAIL
– The linear solver’s initialization function failed.IDA_LSOLVE_FAIL
– The linear solver’s solve function failed in an unrecoverable manner.IDA_BAD_EWT
– Some component of the error weight vector is zero (illegal), either for the input value ofy0
or a corrected value.IDA_FIRST_RES_FAIL
– The user’s residual function returned a recoverable error flag on the first call, butIDACalcIC()
was unable to recover.IDA_RES_FAIL
– The user’s residual function returned a nonrecoverable error flag.IDA_NO_RECOVERY
– The user’s residual function, or the linear solver’s setup or solve function had a recoverable error, butIDACalcIC()
was unable to recover.IDA_CONSTR_FAIL
–IDACalcIC()
was unable to find a solution satisfying the inequality constraints.IDA_LINESEARCH_FAIL
– The linesearch algorithm failed to find a solution with a step larger thansteptol
in weighted RMS norm, and within the allowed number of backtracks.IDA_CONV_FAIL
–IDACalcIC()
failed to get convergence of the Newton iterations.
 Notes:
IDACalcIC()
will correct the values of \(y(t_0)\) and \(\dot{y}(t_0)\) which were specified in the previous call toIDAInit()
orIDAReInit()
. To obtain the corrected values, callIDAGetConsistentIC()
.
6.4.5.8. Rootfinding initialization function
While solving the IVP, IDA has the capability to find the roots of a set of
userdefined functions. To activate the root finding algorithm, call the
following function. This is normally called only once, prior to the first call
to IDASolve()
, but if the rootfinding problem is to be changed during
the solution, IDARootInit()
can also be called prior to a continuation
call to IDASolve()
.

int IDARootInit(void *ida_mem, int nrtfn, IDARootFn g)
The function
IDARootInit
specifies that the roots of a set of functions \(g_i(t,y)\) are to be found while the IVP is being solved. Arguments:
ida_mem
– pointer to the IDA solver objectnrtfn
– is the number of root functionsg
– is the function which defines thenrtfn
functions \(g_i(t,y,\dot{y})\) whose roots are sought. SeeIDARootFn
for more details.
 Return value:
IDA_SUCCESS
– The call was successfulIDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_MEM_FAIL
– A memory allocation failed.IDA_ILL_INPUT
– The functiong
isNULL
, butnrtfn > 0
.
 Notes:
If a new IVP is to be solved with a call to
IDAReInit()
, where the new IVP has no rootfinding problem but the prior one did, then callIDARootInit()
withnrtfn = 0
.
6.4.5.9. IDA solver function
This is the central step in the solution process, the call to perform the
integration of the DAE. The input arguments (itask
) specifies one of two
modes as to where IDA is to return a solution. These modes are modified if
the user has set a stop time (with IDASetStopTime()
) or requested
rootfinding (with IDARootInit()
).

int IDASolve(void *ida_mem, realtype tout, realtype tret, N_Vector yret, N_Vector ypret, int itask)
The function
IDASolve
integrates the DAE over an interval in t. Arguments:
ida_mem
– pointer to the IDA solver object.tout
– the next time at which a computed solution is desired.tret
– the time reached by the solver output.yret
– the computed solution vector y.ypret
– the computed solution vector \(\dot{y}\).itask
– a flag indicating the job of the solver for the next user stepIDA_NORMAL
– the solver will take internal steps until it has reached or just passed the user specifiedtout
parameter. The solver then interpolates in order to return approximate values of \(y(t_{out})\) and \(\dot{y}(t_{out})\).IDA_ONE_STEP
– the solver will just take one internal step and return the solution at the point reached by that step.
 Return value:
IDA_SUCCESS
– The call was successful.IDA_TSTOP_RETURN
–IDASolve()
succeeded by reaching the stop point specified through the optional input functionIDASetStopTime()
.IDA_ROOT_RETURN
–IDASolve()
succeeded and found one or more roots. In this case,tret
is the location of the root. Ifnrtfn
>1, callIDAGetRootInfo()
to see which \(g_i\) were found to have a root.IDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_ILL_INPUT
– One of the inputs toIDASolve()
was illegal, or some other input to the solver was either illegal or missing. The latter category includes the following situations:The tolerances have not been set.
A component of the error weight vector became zero during internal timestepping.
The linear solver initialization function called by the user after calling
IDACreate()
failed to set the linear solverspecificlsolve
field inida_mem
.A root of one of the root functions was found both at a point \(t\) and also very near \(t\).
In any case, the user should see the printed error message for details.
IDA_TOO_MUCH_WORK
– The solver tookmxstep
internal steps but could not reachtout
. The default value formxstep
isMXSTEP_DEFAULT = 500
.IDA_TOO_MUCH_ACC
– The solver could not satisfy the accuracy demanded by the user for some internal step.IDA_ERR_FAIL
– Error test failures occurred too many times (MXNEF = 10
) during one internal time step or occurred with \(h = h_{\text{min}}\).IDA_CONV_FAIL
– Convergence test failures occurred too many times (MXNCF = 10
) during one internal time step or occurred with \(h = h_{\text{min}}\).IDA_LINIT_FAIL
– The linear solver’s initialization function failed.IDA_LSETUP_FAIL
– The linear solver’s setup function failed in an unrecoverable manner.IDA_LSOLVE_FAIL
– The linear solver’s solve function failed in an unrecoverable manner.IDA_CONSTR_FAIL
– The inequality constraints were violated and the solver was unable to recover.IDA_REP_RES_ERR
– The user’s residual function repeatedly returned a recoverable error flag, but the solver was unable to recover.IDA_RES_FAIL
– The user’s residual function returned a nonrecoverable error flag.IDA_RTFUNC_FAIL
– The rootfinding function failed.
 Notes:
The vectors
yret
andypret
can occupy the same space as the initial condition vectorsy0
andyp0
, respectively, that were passed toIDAInit()
.In the
IDA_ONE_STEP
mode,tout
is used on the first call only, and only to get the direction and rough scale of the independent variable.If a stop time is enabled (through a call to
IDASetStopTime()
), thenIDASolve()
returns the solution attstop
. Once the integrator returns at a stop time, any future testing fortstop
is disabled (and can be reenabled only though a new call toIDASetStopTime()
).All failure return values are negative and therefore a test
flag < 0
will trap allIDASolve()
failures.On any error return in which one or more internal steps were taken by
IDASolve()
, the returned values oftret
,yret
, andypret
correspond to the farthest point reached in the integration. On all other error returns, these values are left unchanged from the previousIDASolve()
return.
6.4.5.10. Optional input functions
There are numerous optional input parameters that control the behavior of the IDA solver. IDA provides functions that can be used to change these optional input parameters from their default values. The main inputs are divided in the following categories:
Table 6.1 list the main IDA optional input functions,
Table 6.2 lists the IDALS linear solver interface optional input functions,
Table 6.3 lists the IDANLS nonlinear solver interface optional input functions,
Table 6.4 lists the initial condition calculation optional input functions,
Table 6.5 lists the IDA step size adaptivity optional input functions, and
Table 6.6 lists the rootfinding optional input functions.
These optional inputs are described in detail in the remainder of this section. For the most casual use of IDA, the reader can skip to §6.4.6.
We note that, on an error return, all of the optional input functions also send
an error message to the error handler function. All error return values are
negative, so the test flag < 0
will catch all errors.
The optional input calls can, unless otherwise noted, be executed in any order.
However, if the user’s program calls either IDASetErrFile()
or
IDASetErrHandlerFn()
, then that call should appear first, in order to
take effect for any later error message. Finally, a call to an IDASet***
function can, unless otherwise noted, be made at any time from the user’s
calling program and, if successful, takes effect immediately.
6.4.5.10.1. Main solver optional input functions
Optional input 
Function name 
Default 
Pointer to an error file 


Error handler function 
internal fn. 

User data 
NULL 

Maximum order for BDF method 
5 

Maximum no. of internal steps before \(t_{{\scriptsize out}}\) 
500 

Initial step size 
estimated 

Minimum absolute step size \(h_{\text{min}}\) 
0 

Maximum absolute step size \(h_{\text{max}}\) 
\(\infty\) 

Value of \(t_{stop}\) 
\(\infty\) 

Maximum no. of error test failures 
10 

Suppress alg. vars. from error test 


Variable types (differential/algebraic) 
NULL 

Inequality constraints on solution 
NULL 

int IDASetErrFile(void *ida_mem, FILE *errfp)
The function
IDASetErrFile
specifies the file pointer where all IDA messages should be directed when using the default IDA error handler function. Arguments:
ida_mem
– pointer to the IDA solver object.errfp
– pointer to output file.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The default value for
errfp
isstderr
. Passing a valueNULL
disables all future error message output (except for the case in which the IDA memory pointer isNULL
). This use ofIDASetErrFile()
is strongly discouraged.
Warning
If
IDASetErrFile()
is to be called, it should be called before any other optional input functions, in order to take effect for any later error message.

int IDASetErrHandlerFn(void *ida_mem, IDAErrHandlerFn ehfun, void *eh_data)
The function
IDASetErrHandlerFn
specifies the optional userdefined function to be used in handling error messages. Arguments:
ida_mem
– pointer to the IDA solver object.ehfun
– is the user’s error handler function. SeeIDAErrHandlerFn
for more details.eh_data
– pointer to user data passed toehfun
every time it is called.
 Return value:
IDA_SUCCESS
– The functionehfun
and data pointereh_data
have been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
Error messages indicating that the IDA solver memory is
NULL
will always be directed tostderr
.

int IDASetUserData(void *ida_mem, void *user_data)
The function
IDASetUserData
attaches a userdefined data pointer to the main IDA solver object. Arguments:
ida_mem
– pointer to the IDA solver object.user_data
– pointer to the user data.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
If specified, the pointer to
user_data
is passed to all usersupplied functions that have it as an argument. Otherwise, aNULL
pointer is passed.
Warning
If
user_data
is needed in user linear solver or preconditioner functions, the call toIDASetUserData()
must be made before the call to specify the linear solver.

int IDASetMaxOrd(void *ida_mem, int maxord)
The function
IDASetMaxOrd
specifies the maximum order of the linear multistep method. Arguments:
ida_mem
– pointer to the IDA solver object.maxord
– value of the maximum method order. This must be positive.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– The input valuemaxord
is \(\leq\) 0 , or larger than the max order value whenIDAInit()
was called.
 Notes:
The default value is 5. If the input value exceeds 5, the value 5 will be used. If called before
IDAInit()
,maxord
limits the memory requirements for the internal IDA memory block and its value cannot be increased past the value set whenIDAInit()
was called.

int IDASetMaxNumSteps(void *ida_mem, long int mxsteps)
The function
IDASetMaxNumSteps
specifies the maximum number of steps to be taken by the solver in its attempt to reach the next output time. Arguments:
ida_mem
– pointer to the IDA solver object.mxsteps
– maximum allowed number of steps.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
Passing
mxsteps
= 0 results in IDA using the default value (500). Passingmxsteps
< 0 disables the test (not recommended).

int IDASetInitStep(void *ida_mem, realtype hin)
The function
IDASetInitStep
specifies the initial step size. Arguments:
ida_mem
– pointer to the IDA solver object.hin
– value of the initial step size to be attempted. Pass 0.0 to have IDA use the default value.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
By default, IDA estimates the initial step as the solution of \(\h \dot{y} \_{{\scriptsize WRMS}} = 1/2\), with an added restriction that \(h \leq .001t_{\text{out}}  t_0\).

int IDASetMinStep(void *ida_mem, realtype hmin)
The function
IDASetMinStep
specifies the minimum absolute value of the step size.Pass
hmin = 0
to obtain the default value of 0. Arguments:
ida_mem
– pointer to the IDA solver object.hmin
– minimum absolute value of the step size.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
–hmin
is negative.
New in version 6.2.0.

int IDASetMaxStep(void *ida_mem, realtype hmax)
The function
IDASetMaxStep
specifies the maximum absolute value of the step size. Arguments:
ida_mem
– pointer to the IDA solver object.hmax
– maximum absolute value of the step size.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– Eitherhmax
is not positive or it is smaller than the minimum allowable step.
 Notes:
Pass
hmax = 0
to obtain the default value \(\infty\).

int IDASetStopTime(void *ida_mem, realtype tstop)
The function
IDASetStopTime
specifies the value of the independent variable \(t\) past which the solution is not to proceed. Arguments:
ida_mem
– pointer to the IDA solver object.tstop
– value of the independent variable past which the solution should not proceed.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– The value oftstop
is not beyond the current \(t\) value, \(t_n\).
 Notes:
The default, if this routine is not called, is that no stop time is imposed. Once the integrator returns at a stop time, any future testing for
tstop
is disabled (and can be reenabled only though a new call toIDASetStopTime()
).

int IDASetMaxErrTestFails(void *ida_mem, int maxnef)
The function
IDASetMaxErrTestFails
specifies the maximum number of error test failures in attempting one step. Arguments:
ida_mem
– pointer to the IDA solver object.maxnef
– maximum number of error test failures allowed on one step (>0).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The default value is 10.

int IDASetSuppressAlg(void *ida_mem, booleantype suppressalg)
The function
IDASetSuppressAlg
indicates whether or not to suppress algebraic variables in the local error test. Arguments:
ida_mem
– pointer to the IDA solver object.suppressalg
– indicates whether to suppress (SUNTRUE
) or include (SUNFALSE
) the algebraic variables in the local error test.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The default value is
SUNFALSE
. Ifsuppressalg = SUNTRUE
is selected, then theid
vector must be set (throughIDASetId()
) to specify the algebraic components. In general, the use of this option (withsuppressalg = SUNTRUE
) is discouraged when solving DAE systems of index 1, whereas it is generally encouraged for systems of index 2 or more. See pp. 146147 of [19] for more on this issue.

int IDASetId(void *ida_mem, N_Vector id)
The function
IDASetId
specifies algebraic/differential components in the \(y\) vector. Arguments:
ida_mem
– pointer to the IDA solver object.id
– a vector of values identifying the components of \(y\) as differential or algebraic variables. A value of 1.0 indicates a differential variable, while 0.0 indicates an algebraic variable.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The vector
id
is required if the algebraic variables are to be suppressed from the local error test (seeIDASetSuppressAlg()
) or ifIDACalcIC()
is to be called withicopt
=IDA_YA_YDP_INIT
.

int IDASetConstraints(void *ida_mem, N_Vector constraints)
The function
IDASetConstraints
specifies a vector defining inequality constraints for each component of the solution vector \(y\). Arguments:
ida_mem
– pointer to the IDA solver object.constraints
– vector of constraint flags.If
constraints[i] = 0
, no constraint is imposed on \(y_i\).If
constraints[i] = 1
, \(y_i\) will be constrained to be \(y_i \ge 0.0\).If
constraints[i] = 1
, \(y_i\) will be constrained to be \(y_i \le 0.0\).If
constraints[i] = 2
, \(y_i\) will be constrained to be \(y_i > 0.0\).If
constraints[i] = 2
, \(y_i\) will be constrained to be \(y_i < 0.0\).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– The constraints vector contains illegal values.
 Notes:
The presence of a non
NULL
constraints vector that is not 0.0 in all components will cause constraint checking to be performed. However, a call with 0.0 in all components of constraints vector will result in an illegal input return. ANULL
input will disable constraint checking.
6.4.5.10.2. Linear solver interface optional input functions
Optional input 
Function name 
Default 
Jacobian function 
DQ 

Set parameter determining if a \(c_j\) change requires a linear solver setup call 
0.25 

Enable or disable linear solution scaling 
on 

Jacobiantimesvector function 
NULL, DQ 

Preconditioner functions 
NULL, NULL 

Ratio between linear and nonlinear tolerances 
0.05 

Increment factor used in DQ \(Jv\) approx. 
1.0 

Jacobiantimesvector DQ Res function 
NULL 

Newton linear solve tolerance conversion factor 
vector length 
The mathematical explanation of the linear solver methods available to IDA is provided in §6.2.2. We group the usercallable routines into four categories: general routines concerning the overall IDALS linear solver interface, optional inputs for matrixbased linear solvers, optional inputs for matrixfree linear solvers, and optional inputs for iterative linear solvers. We note that the matrixbased and matrixfree groups are mutually exclusive, whereas the “iterative” tag can apply to either case.
When using matrixbased linear solver modules, the IDALS solver interface needs
a function to compute an approximation to the Jacobian matrix
\(J(t,y,\dot{y})\). This function must be of type IDALsJacFn
. The
user can supply a Jacobian function or, if using a
SUNMATRIX_DENSE or SUNMATRIX_BAND matrix
\(J\), can use the default internal difference quotient approximation that
comes with the IDALS interface. To specify a usersupplied Jacobian function
jac
, IDALS provides the function IDASetJacFn()
. The IDALS interface
passes the pointer user_data
to the Jacobian function. This allows the user
to create an arbitrary structure with relevant problem data and access it during
the execution of the usersupplied Jacobian function, without using global data
in the program. The pointer user_data
may be specified through
IDASetUserData()
.

int IDASetJacFn(void *ida_mem, IDALsJacFn jac)
The function
IDASetJacFn
specifies the Jacobian approximation function to be used for a matrixbased solver within the IDALS interface. Arguments:
ida_mem
– pointer to the IDA solver object.jac
– userdefined Jacobian approximation function. SeeIDALsJacFn
for more details.
 Return value:
IDALS_SUCCESS
– The optional value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver interface has not been initialized.
 Notes:
This function must be called after the IDALS linear solver interface has been initialized through a call to
IDASetLinearSolver()
. By default, IDALS uses an internal difference quotient function for the SUNMATRIX_DENSE and SUNMATRIX_BAND modules. IfNULL
is passed tojac
, this default function is used. An error will occur if nojac
is supplied when using other matrix types.
Warning
The previous routine
IDADlsSetJacFn
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
When using a matrixbased linear solver the matrix information will be updated
infrequently to reduce matrix construction and, with direct solvers,
factorization costs. As a result the value of \(\alpha\) may not be current
and a scaling factor is applied to the solution of the linear system to account
for the lagged value of \(\alpha\). See §11.5.1 for
more details. The function IDASetLinearSolutionScaling()
can be used to
disable this scaling when necessary, e.g., when providing a custom linear solver
that updates the matrix using the current \(\alpha\) as part of the solve.

int IDASetDeltaCjLSetup(void *ida_mem, realtype dcj)
The function
IDASetDeltaCjLSetup
specifies the parameter that determines the bounds on a change in \(c_j\) that require a linear solver setup call. Ifcj_current / cj_previous < (1  dcj) / (1 + dcj)
orcj_current / cj_previous > (1 + dcj) / (1  dcj)
, the linear solver setup function is called.If
dcj
is \(< 0\) or \(\geq 1\) then the default value (0.25) is used. Arguments:
ida_mem
– pointer to the IDA memory block.dcj
– the \(c_j\) change threshold.
 Return value:
IDA_SUCCESS
– The flag value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetLinearSolutionScaling(void *ida_mem, booleantype onoff)
The function
IDASetLinearSolutionScaling
enables or disables scaling the linear system solution to account for a change in \(\alpha\) in the linear system. For more details see §11.5.1. Arguments:
ida_mem
– pointer to the IDA solver object.onoff
– flag to enable (SUNTRUE
) or disable (SUNFALSE
) scaling
 Return value:
IDALS_SUCCESS
– The flag value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver interface has not been initialized.IDALS_ILL_INPUT
– The attached linear solver is not matrixbased.
 Notes:
This function must be called after the IDALS linear solver interface has been initialized through a call to
IDASetLinearSolver()
. By default scaling is enabled with matrixbased linear solvers.
When using matrixfree linear solver modules, the IDALS solver interface requires a function to compute an approximation to the product between the Jacobian matrix \(J(t,y,\dot{y})\) and a vector \(v\). The user can supply a Jacobiantimesvector approximation function, or use the default internal difference quotient function that comes with the IDALS solver interface.
A userdefined Jacobianvector product function must be of type
IDALsJacTimesVecFn
and can be specified through a call to
IDASetJacTimes()
. The evaluation and processing of any Jacobianrelated
data needed by the user’s Jacobianvector product function may be done in the
optional usersupplied function jtsetup
(see
§6.4.6.7 for specification details). The
pointer user_data
received through IDASetUserData()
(or a pointer to
NULL
if user_data
was not specified) is passed to the Jacobianvector
product setup and product functions, jtsetup
and jtimes
, each time they
are called. This allows the user to create an arbitrary structure with relevant
problem data and access it during the execution of the usersupplied functions
without using global data in the program.

int IDASetJacTimes(void *ida_mem, IDALsJacTimesSetupFn jsetup, IDALsJacTimesVecFn jtimes)
The function
IDASetJacTimes
specifies the Jacobianvector product setup and product functions. Arguments:
ida_mem
– pointer to the IDA solver object.jtsetup
– userdefined function to set up the Jacobianvector product. SeeIDALsJacTimesSetupFn
for more details. PassNULL
if no setup is necessary.jtimes
– userdefined Jacobianvector product function. SeeIDALsJacTimesVecFn
for more details.
 Return value:
IDALS_SUCCESS
– The optional value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.IDALS_SUNLS_FAIL
– An error occurred when setting up the system matrixtimesvector routines in theSUNLinearSolver
object used by the IDALS interface.
 Notes:
The default is to use an internal finite difference quotient for
jtimes
and to omitjtsetup
. IfNULL
is passed tojtimes
, these defaults are used. A user may specify nonNULL
jtimes
andNULL
jtsetup
inputs. This function must be called after the IDALS linear solver interface has been initialized through a call toIDASetLinearSolver()
.
Warning
The previous routine
IDASpilsSetJacTimes
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
When using the default differencequotient approximation to the Jacobianvector
product, the user may specify the factor to use in setting increments for the
finitedifference approximation, via a call to IDASetIncrementFactor()
.

int IDASetIncrementFactor(void *ida_mem, realtype dqincfac)
The function
IDASetIncrementFactor
specifies the increment factor to be used in the differencequotient approximation to the product \(Jv\). Specifically, \(Jv\) is approximated via the formula\[Jv = \frac{1}\sigma\left[F(t,\tilde{y},\tilde{\dot{y}})  F(t,y,\dot{y})\right],\]where \(\tilde{y} = y + \sigma v\), \(\tilde{\dot{y}} = \dot{y} + c_j \sigma v\), \(c_j\) is a BDF parameter proportional to the step size, \(\sigma = \mathtt{dqincfac} \sqrt{N}\), and \(N\) is the number of equations in the DAE system.
 Arguments:
ida_mem
– pointer to the IDA solver object.dqincfac
– userspecified increment factor positive.
 Return value:
IDALS_SUCCESS
– The optional value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.IDALS_ILL_INPUT
– The specified value ofdqincfac
is \(\le 0\).
 Notes:
The default value is 1.0. This function must be called after the IDALS linear solver interface has been initialized through a call to
IDASetLinearSolver()
.
Warning
The previous routine
IDASpilsSetIncrementFactor()
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
Additionally, when using the internal difference quotient, the user may also
optionally supply an alternative residual function for use in the
Jacobianvector product approximation by calling
IDASetJacTimesResFn()
. The alternative residual function should compute
a suitable (and differentiable) approximation to the residual function provided
to IDAInit()
. For example, as done in [44] for an
ODE in explicit form, the alternative function may use lagged values when
evaluating a nonlinearity to avoid differencing a potentially nondifferentiable
factor.

int IDASetJacTimesResFn(void *ida_mem, IDAResFn jtimesResFn)
The function
IDASetJacTimesResFn
specifies an alternative DAE residual function for use in the internal Jacobianvector product difference quotient approximation. Arguments:
ida_mem
– pointer to the IDA solver object.jtimesResFn
– is the function which computes the alternative DAE residual function to use in Jacobianvector product difference quotient approximations. SeeIDAResFn
for more details.
 Return value:
IDALS_SUCCESS
– The optional value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.IDALS_ILL_INPUT
– The internal difference quotient approximation is disabled.
 Notes:
The default is to use the residual function provided to
IDAInit()
in the internal difference quotient. If the input resudual function isNULL
, the default is used. This function must be called after the IDALS linear solver interface has been initialized through a call toIDASetLinearSolver()
.
When using an iterative linear solver, the user may supply a preconditioning
operator to aid in solution of the system. This operator consists of two
usersupplied functions, psetup
and psolve
, that are supplied to IDA
using the function IDASetPreconditioner()
. The psetup
function
supplied to this routine should handle evaluation and preprocessing of any
Jacobian data needed by the user’s preconditioner solve function,
psolve
. Both of these functions are fully specified in
§6.4.6.8 and
§6.4.6.9). The user data pointer received
through IDASetUserData()
(or NULL
if a user data pointer was not
specified) is passed to the psetup
and psolve
functions. This allows the
user to create an arbitrary structure with relevant problem data and access it
during the execution of the usersupplied preconditioner functions without using
global data in the program.

int IDASetPreconditioner(void *ida_mem, IDALsPrecSetupFn psetup, IDALsPrecSolveFn psolve)
The function
IDASetPreconditioner
specifies the preconditioner setup and solve functions. Arguments:
ida_mem
– pointer to the IDA solver object.psetup
– userdefined function to set up the preconditioner. SeeIDALsPrecSetupFn
for more details. PassNULL
if no setup is necessary.psolve
– userdefined preconditioner solve function. SeeIDALsPrecSolveFn
for more details.
 Return value:
IDALS_SUCCESS
– The optional values have been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.IDALS_SUNLS_FAIL
– An error occurred when setting up preconditioning in theSUNLinearSolver
object used by the IDALS interface.
 Notes:
The default is
NULL
for both arguments (i.e., no preconditioning). This function must be called after the IDALS linear solver interface has been initialized through a call toIDASetLinearSolver()
.
Warning
The previous routine
IDASpilsSetPreconditioner
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
Also, as described in §6.2.2, the IDALS interface requires that iterative linear solvers stop when the norm of the preconditioned residual satisfies
where \(\epsilon\) is the nonlinear solver tolerance, and the default
\(\epsilon_L = 0.05\); this value may be modified by the user through the
IDASetEpsLin()
function.

int IDASetEpsLin(void *ida_mem, realtype eplifac)
The function
IDASetEpsLin
specifies the factor by which the Krylov linear solver’s convergence test constant is reduced from the nonlinear iteration test constant. Arguments:
ida_mem
– pointer to the IDA solver object.eplifac
– linear convergence safety factor \(\geq 0.0\).
 Return value:
IDALS_SUCCESS
– The optional value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.IDALS_ILL_INPUT
– The factoreplifac
is negative.
 Notes:
The default value is \(0.05\). This function must be called after the IDALS linear solver interface has been initialized through a call to
IDASetLinearSolver()
. Ifeplifac
\(= 0.0\) is passed, the default value is used.
Warning
The previous routine
IDASpilsSetEpsLin
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDASetLSNormFactor(void *ida_mem, realtype nrmfac)
The function
IDASetLSNormFactor
specifies the factor to use when converting from the integrator tolerance (WRMS norm) to the linear solver tolerance (L2 norm) for Newton linear system solves e.g.,tol_L2 = fac * tol_WRMS
. Arguments:
ida_mem
– pointer to the IDA solver object.nrmfac
– the norm conversion factor.If
nrmfac > 0
, the provided value is used.If
nrmfac = 0
then the conversion factor is computed using the vector length i.e.,nrmfac = N_VGetLength(y)
(default).If
nrmfac < 0
then the conversion factor is computed using the vector dot productnrmfac = N_VDotProd(v,v)
where all the entries ofv
are one.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
This function must be called after the IDALS linear solver interface has been initialized through a call to
IDASetLinearSolver()
. Prior to the introduction ofN_VGetLength()
in SUNDIALS v5.0.0 (IDA v5.0.0) the value ofnrmfac
was computed usingN_VDotProd()
i.e., thenrmfac < 0
case.
6.4.5.10.3. Nonlinear solver interface optional input functions
Optional input 
Function name 
Default 
Maximum no. of nonlinear iterations 
4 

Maximum no. of convergence failures 
10 

Coeff. in the nonlinear convergence test 
0.33 

Residual function for nonlinear system evaluations 

The following functions can be called to set optional inputs controlling the nonlinear solver.

int IDASetMaxNonlinIters(void *ida_mem, int maxcor)
The function
IDASetMaxNonlinIters
specifies the maximum number of nonlinear solver iterations in one solve attempt. Arguments:
ida_mem
– pointer to the IDA solver object.maxcor
– maximum number of nonlinear solver iterations allowed in one solve attempt (>0).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_MEM_FAIL
– TheSUNNonlinearSolver
object isNULL
.
 Notes:
The default value is 4.

int IDASetMaxConvFails(void *ida_mem, int maxncf)
The function
IDASetMaxConvFails
specifies the maximum number of nonlinear solver convergence failures in one step. Arguments:
ida_mem
– pointer to the IDA solver object.maxncf
– maximum number of allowable nonlinear solver convergence failures in one step (>0).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The default value is 10.

int IDASetNonlinConvCoef(void *ida_mem, realtype nlscoef)
The function
IDASetNonlinConvCoef
specifies the safety factor in the nonlinear convergence test; see (6.7). Arguments:
ida_mem
– pointer to the IDA solver object.nlscoef
– coefficient in nonlinear convergence test (>0.0).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– The value ofnlscoef
is \(\leq 0.0\).
 Notes:
The default value is 0.33.

int IDASetNlsResFn(void *ida_mem, IDAResFn res)
The function
IDASetNlsResFn
specifies an alternative residual function for use in nonlinear system function evaluations. Arguments:
ida_mem
– pointer to the IDA solver object.res
– the alternative function which computes the DAE residual function \(F(t, y, \dot{y})\). SeeIDAResFn
for more details.
 Return value:
IDA_SUCCESS
– The optional function has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The default is to use the residual function provided to
IDAInit()
in nonlinear system function evaluations. If the input residual function isNULL
, the default is used.When using a nondefault nonlinear solver, this function must be called after
IDASetNonlinearSolver()
.
6.4.5.10.4. Initial condition calculation optional input functions
Optional input 
Function name 
Default 
Coeff. in the nonlinear convergence test 
0.0033 

Maximum no. of steps 
5 

Maximum no. of Jacobian/precond. evals. 
4 

Maximum no. of Newton iterations 
10 

Max. linesearch backtracks per Newton iter. 
100 

Turn off linesearch 


Lower bound on Newton step 
uround\(^{2/3}\) 
The following functions can be called just prior to calling IDACalcIC()
to set optional inputs controlling the initial condition calculation.

int IDASetNonlinConvCoefIC(void *ida_mem, realtype epiccon)
The function
IDASetNonlinConvCoefIC
specifies the positive constant in the Newton iteration convergence test within the initial condition calculation. Arguments:
ida_mem
– pointer to the IDA solver object.epiccon
– coefficient in the Newton convergence test \((>0)\).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– Theepiccon
factor is \(\leq 0.0\).
 Notes:
The default value is \(0.01 \cdot 0.33\). This test uses a weighted RMS norm (with weights defined by the tolerances). For new initial value vectors \(y\) and \(\dot{y}\) to be accepted, the norm of \(J^{1}F(t_0, y, \dot{y})\) must be \(\leq \mathtt{epiccon}\), where \(J\) is the system Jacobian.

int IDASetMaxNumStepsIC(void *ida_mem, int maxnh)
The function
IDASetMaxNumStepsIC
specifies the maximum number of steps allowed whenicopt = IDA_YA_YDP_INIT
inIDACalcIC()
, where \(h\) appears in the system Jacobian, \(J = \dfrac{\partial F}{\partial y} + \left(\dfrac1h\right)\dfrac{\partial F}{\partial \dot{y}}\). Arguments:
ida_mem
– pointer to the IDA solver object.maxnh
– maximum allowed number of values for \(h\).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
–maxnh
is nonpositive.
 Notes:
The default value is \(5\).

int IDASetMaxNumJacsIC(void *ida_mem, int maxnj)
The function
IDASetMaxNumJacsIC
specifies the maximum number of the approximate Jacobian or preconditioner evaluations allowed when the Newton iteration appears to be slowly converging. Arguments:
ida_mem
– pointer to the IDA solver object.maxnj
– maximum allowed number of Jacobian or preconditioner evaluations.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
–maxnj
is nonpositive.
 Notes:
The default value is \(4\).

int IDASetMaxNumItersIC(void *ida_mem, int maxnit)
The function
IDASetMaxNumItersIC
specifies the maximum number of Newton iterations allowed in any one attempt to solve the initial conditions calculation problem. Arguments:
ida_mem
– pointer to the IDA solver object.maxnit
– maximum number of Newton iterations.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
–maxnit
is nonpositive.
 Notes:
The default value is \(10\).

int IDASetMaxBacksIC(void *ida_mem, int maxbacks)
The function
IDASetMaxBacksIC
specifies the maximum number of linesearch backtracks allowed in any Newton iteration, when solving the initial conditions calculation problem. Arguments:
ida_mem
– pointer to the IDA solver object.maxbacks
– maximum number of linesearch backtracks per Newton step.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
–maxbacks
is nonpositive.
 Notes:
The default value is \(100\).

int IDASetLineSearchOffIC(void *ida_mem, booleantype lsoff)
The function
IDASetLineSearchOffIC
specifies whether to turn on or off the linesearch algorithm. Arguments:
ida_mem
– pointer to the IDA solver object.lsoff
– a flag to turn off (SUNTRUE
) or keep (SUNFALSE
) the linesearch algorithm.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
Notes:
The default value is
SUNFALSE
.

int IDASetStepToleranceIC(void *ida_mem, int steptol)
The function
IDASetStepToleranceIC
specifies a positive lower bound on the Newton step. Arguments:
ida_mem
– pointer to the IDA solver object.steptol
– Minimum allowed WRMSnorm of the Newton step \((> 0.0)\).
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– Thesteptol
tolerance is \(\leq 0.0\).
 Notes:
The default value is \((\text{unit roundoff})^{2/3}\).
6.4.5.10.5. Time step adaptivity optional input functions
Optional input 
Function name 
Default 

Fixed step size bounds \(\eta_{\mathrm{min\_fx}}\) and \(\eta_{\mathrm{max\_fx}}\) 
1.0 and 2.0 

Maximum step size growth factor \(\eta_{\mathrm{max}}\) 
2.0 

Minimum step size reduction factor \(\eta_{\mathrm{min}}\) 
0.5 

Maximum step size reduction factor \(\eta_{\mathrm{low}}\) 
0.9 

Minimum step size reduction factor after an error test failure \(\eta_{\mathrm{min\_ef}}\) 
0.25 

Step size reduction factor after a nonlinear solver convergence failure \(\eta_{\mathrm{cf}}\) 
0.25 
The following functions can be called to set optional inputs to control the step size adaptivity.
Note
The default values for the step size adaptivity tuning parameters have a long history of success and changing the values is generally discouraged. However, users that wish to experiment with alternative values should be careful to make changes gradually and with testing to determine their effectiveness.

int IDASetEtaFixedStepBounds(void *ida_mem, realtype eta_min_fx, realtype eta_max_fx)
The function
IDASetEtaFixedStepBounds
specifies the bounds \(\eta_{\mathrm{min\_fx}}\) and \(\eta_{\mathrm{max\_fx}}\). If step size change factor \(\eta\) satisfies \(\eta_{\mathrm{min\_fx}} < \eta < \eta_{\mathrm{max\_fx}}\) the current step size is retained.The default values are \(\eta_{\mathrm{fxmin}} = 1\) and \(\eta_{\mathrm{fxmax}} = 2\).
eta_fxmin
should satisfy \(0 < \eta_{\mathrm{fxmin}} \leq 1\), otherwise the default value is used.eta_fxmax
should satisfy \(\eta_{\mathrm{fxmin}} \geq 1\), otherwise the default value is used. Arguments:
ida_mem
– pointer to the IDA solver object.eta_min_fx
– value of the fixed step size lower bound.eta_max_fx
– value of the fixed step size upper bound.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetEtaMax(void *ida_mem, realtype eta_max)
The function
IDASetEtaMax
specifies the maximum step size growth factor \(\eta_{\mathrm{max}}\).The default value is \(\eta_{\mathrm{max}} = 2\).
 Arguments:
ida_mem
– pointer to the IDA solver object.eta_max
– maximum step size growth factor.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetEtaMin(void *ida_mem, realtype eta_min)
The function
IDASetEtaMin
specifies the minimum step size reduction factor \(\eta_{\mathrm{min}}\).The default value is \(\eta_{\mathrm{min}} = 0.5\).
eta_min
should satisfy \(0 < \eta_{\mathrm{min}} < 1\), otherwise the default value is used. Arguments:
ida_mem
– pointer to the IDA solver object.eta_min
– value of the minimum step size reduction factor.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetEtaLow(void *ida_mem, realtype eta_low)
The function
IDASetEtaLow
specifies the maximum step size reduction factor \(\eta_{\mathrm{low}}\).The default value is \(\eta_{\mathrm{low}} = 0.9\).
eta_low
should satisfy \(0 < \eta_{\mathrm{low}} \leq 1\), otherwise the default value is used. Arguments:
ida_mem
– pointer to the IDA solver object.eta_low
– value of the maximum step size reduction factor.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetEtaMinErrFail(void *ida_mem, realtype eta_min_ef)
The function
IDASetEtaMinErrFail
specifies the minimum step size reduction factor \(\eta_{\mathrm{min\_ef}}\) after an error test failure.The default value is \(\eta_{\mathrm{min\_ef}} = 0.25\).
If
eta_min_ef
is \(\leq 0\) or \(\geq 1\), the default value is used. Arguments:
ida_mem
– pointer to the IDA solver object.eta_low
– value of the minimum step size reduction factor.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.

int IDASetEtaConvFail(void *ida_mem, realtype eta_cf)
The function
IDASetEtaConvFail
specifies the step size reduction factor \(\eta_{\mathrm{cf}}\) after a nonlinear solver convergence failure.The default value is \(\eta_{\mathrm{cf}} = 0.25\).
If
eta_cf
is \(\leq 0\) or \(\geq 1\), the default value is used. Arguments:
ida_mem
– pointer to the IDA solver object.eta_low
– value of the step size reduction factor.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.2.0.
6.4.5.10.6. Rootfinding optional input functions
Optional input 
Function name 
Default 
Direction of zerocrossing 
both 

Disable rootfinding warnings 
none 
The following functions can be called to set optional inputs to control the rootfinding algorithm.

int IDASetRootDirection(void *ida_mem, int *rootdir)
The function
IDASetRootDirection
specifies the direction of zerocrossings to be located and returned to the user. Arguments:
ida_mem
– pointer to the IDA solver object.rootdir
– state array of lengthnrtfn
, the number of root functions \(g_i\) , as specified in the call to the functionIDARootInit()
.A value of \(0\) for
rootdir[i]
indicates that crossing in either direction should be reported for \(g_i\).A value of \(+1\) or \(1\) for
rootdir[i]
indicates that the solver should report only zerocrossings where \(g_i\) is increasing or decreasing, respectively.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– rootfinding has not been activated through a call toIDARootInit()
.
 Notes:
The default behavior is to locate both zerocrossing directions.

int IDASetNoInactiveRootWarn(void *ida_mem)
The function
IDASetNoInactiveRootWarn
disables issuing a warning if some root function appears to be identically zero at the beginning of the integration. Arguments:
ida_mem
– pointer to the IDA solver object.
 Return value:
IDA_SUCCESS
– The optional value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
IDA will not report the initial conditions as a possible zerocrossing (assuming that one or more components \(g_i\) are zero at the initial time). However, if it appears that some \(g_i\) is identically zero at the initial time (i.e., \(g_i\) is zero at the initial time and after the first step), IDA will issue a warning which can be disabled with this optional input function.
6.4.5.11. Interpolated output function
An optional function IDAGetDky()
is available to obtain additional
output values. This function must be called after a successful return from
IDASolve()
and provides interpolated values of \(y\) or its
derivatives of order up to the last internal order used for any value of
\(t\) in the last internal step taken by IDA.

int IDAGetDky(void *ida_mem, realtype t, int k, N_Vector dky)
The function
IDAGetDky
computes the interpolated values of the \(k^{th}\) derivative of \(y\) for any value of \(t\) in the last internal step taken by IDA. The value of \(k\) must be nonnegative and smaller than the last internal order used. A value of \(0\) for \(k\) means that the \(y\) is interpolated. The value of \(t\) must satisfy \(t_n  h_u \le t \le t_n\), where \(t_n\) denotes the current internal time reached, and \(h_u\) is the last internal step size used successfully. Arguments:
ida_mem
– pointer to the IDA solver object.t
– time at which to interpolate.k
– integer specifying the order of the derivative of \(y\) wanted.dky
– vector containing the interpolated \(k^{th}\) derivative of \(y(t)\).
 Return value:
IDA_SUCCESS
–IDAGetDky
succeeded.IDA_MEM_NULL
– Theida_mem
argument wasNULL
.IDA_BAD_T
–t
is not in the interval \([t_n  h_u , t_n]\).IDA_BAD_K
–k
is not one of \({0, 1, \ldots, k_{\text{last}}}\).IDA_BAD_DKY
–dky
isNULL
.
 Notes:
It is only legal to call the function
IDAGetDky()
after a successful return fromIDASolve()
. FunctionsIDAGetCurrentTime()
,IDAGetLastStep()
andIDAGetLastOrder()
can be used to access \(t_n\), \(h_u\), and \(k_{\text{last}}\).
6.4.5.12. Optional output functions
IDA provides an extensive list of functions that can be used to obtain solver performance information. Table 6.7 lists all optional output functions in IDA, which are then described in detail in the remainder of this section.
Some of the optional outputs, especially the various counters, can be very
useful in determining how successful the IDA solver is in doing its job. For
example, the counters nsteps
and nrevals
provide a rough measure of the
overall cost of a given run, and can be compared among runs with differing input
options to suggest which set of options is most efficient. The ratio
nniters/nsteps
measures the performance of the nonlinear solver in solving
the nonlinear systems at each time step; typical values for this range from 1.1
to 1.8. The ratio njevals/nniters
(in the case of a matrixbased linear
solver), and the ratio npevals/nniters
(in the case of an iterative linear
solver) measure the overall degree of nonlinearity in these systems, and also
the quality of the approximate Jacobian or preconditioner being used. Thus, for
example, njevals/nniters
can indicate if a usersupplied Jacobian is
inaccurate, if this ratio is larger than for the case of the corresponding
internal Jacobian. The ratio nliters/nniters
measures the performance of the
Krylov iterative linear solver, and thus (indirectly) the quality of the
preconditioner.
Optional output 
Function name 

Size of IDA real and integer workspace 

Cumulative number of internal steps 

No. of calls to residual function 

No. of calls to linear solver setup function 

No. of local error test failures that have occurred 

No. of failed steps due to a nonlinear solver failure 

Order used during the last step 

Order to be attempted on the next step 

Actual initial step size used 

Step size used for the last step 

Step size to be attempted on the next step 

Current internal time reached by the solver 

Suggested factor for tolerance scaling 

Error weight vector for state variables 

Estimated local errors 

All IDA integrator statistics 

No. of nonlinear solver iterations 

No. of nonlinear convergence failures 

IDA nonlinear solver statistics 

User data pointer 

Array showing roots found 

No. of calls to user root function 

Print all statistics 

Name of constant associated with a return flag 

Number of backtrack operations 

Corrected initial conditions 

Size of real and integer workspace 

No. of Jacobian evaluations 

No. of residual calls for finite diff. Jacobianvector evals. 

No. of linear iterations 

No. of linear convergence failures 

No. of preconditioner evaluations 

No. of preconditioner solves 

No. of Jacobianvector setup evaluations 

No. of Jacobianvector product evaluations 

Last return from a linear solver function 

Name of constant associated with a return flag 
6.4.5.12.1. Main solver optional output functions
IDA provides several usercallable functions that can be used to obtain different quantities that may be of interest to the user, such as solver workspace requirements, solver performance statistics, as well as additional data from the IDA solver object (a suggested tolerance scaling factor, the error weight vector, and the vector of estimated local errors). Also provided are functions to extract statistics related to the performance of the nonlinear solver being used. As a convenience, additional extraction functions provide the optional outputs in groups. These optional output functions are described next.

int IDAGetWorkSpace(void *ida_mem, long int lenrw, long int leniw)
The function
IDAGetWorkSpace
returns the IDA real and integer workspace sizes. Arguments:
ida_mem
– pointer to the IDA solver object.lenrw
– number of real values in the IDA workspace.leniw
– number of integer values in the IDA workspace.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
In terms of the problem size \(N\), the maximum method order
maxord
, and the number of root functionsnrtfn
(see §6.4.5.8), the actual size of the real workspace, inrealtype
words, is given by the following:base value: \(\mathtt{lenrw} = 55 + (m + 6) * N_r + 3 * \mathtt{nrtfn}\);
with
IDASVtolerances()
: \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);with constraint checking (see
IDASetConstraints()
): \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);with
id
specified (seeIDASetId()
): \(\mathtt{lenrw} = \mathtt{lenrw} + N_r\);
where \(m = \max(\mathtt{maxord}, 3)\), and \(N_r\) is the number of real words in one
N_Vector
\((\approx N)\).The size of the integer workspace (without distinction between
int
andlong int
words) is given by:base value: \(\mathtt{leniw} = 38 + (m + 6) * N_i + \mathtt{nrtfn}\);
with
IDASVtolerances()
: \(\mathtt{leniw} = \mathtt{leniw} + N_i\);with constraint checking: \(\mathtt{lenrw} = \mathtt{lenrw} + N_i\);
with
id
specified (seeIDASetId()
): \(\mathtt{lenrw} = \mathtt{lenrw} + N_i\);
where \(N_i\) is the number of integer words in one
N_Vector
(= 1 for the serialN_Vector
and2 * npes
for the parallelN_Vector
onnpes
processors). For the default value ofmaxord
, with no rootfinding, noid
, no constraints, and with no call toIDASVtolerances()
, these lengths are given roughly by \(\mathtt{lenrw} = 55 + 11 * N\) and \(\mathtt{leniw} = 49\).

int IDAGetNumSteps(void *ida_mem, long int *nsteps)
The function
IDAGetNumSteps
returns the cumulative number of internal steps taken by the solver (total so far). Arguments:
ida_mem
– pointer to the IDA solver object.nsteps
– number of steps taken by IDA.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetNumResEvals(void *ida_mem, long int *nrevals)
The function
IDAGetNumResEvals
returns the number of calls to the user’s residual evaluation function. Arguments:
ida_mem
– pointer to the IDA solver object.nrevals
– number of calls to the user’sres
function.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
The
nrevals
value returned byIDAGetNumResEvals()
does not account for calls made tores
from a linear solver or preconditioner module.

int IDAGetNumLinSolvSetups(void *ida_mem, long int *nlinsetups)
The function
IDAGetNumLinSolvSetups
returns the cumulative number of calls made to the linear solver’s setup function (total so far). Arguments:
ida_mem
– pointer to the IDA solver object.nlinsetups
– number of calls made to the linear solver setup function.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetNumErrTestFails(void *ida_mem, long int *netfails)
The function
IDAGetNumErrTestFails
returns the cumulative number of local error test failures that have occurred (total so far). Arguments:
ida_mem
– pointer to the IDA solver object.netfails
– number of error test failures.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetNumStepSolveFails(void *ida_mem, long int *ncnf)
Returns the number of failed steps due to a nonlinear solver failure.
 Arguments:
ida_mem
– pointer to the IDA solver object.ncnf
– number of step failures.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetLastOrder(void *ida_mem, int *klast)
The function
IDAGetLastOrder
returns the integration method order used during the last internal step. Arguments:
ida_mem
– pointer to the IDA solver object.klast
– method order used on the last internal step.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetCurrentOrder(void *ida_mem, int *kcur)
The function
IDAGetCurrentOrder
returns the integration method order to be used on the next internal step. Arguments:
ida_mem
– pointer to the IDA solver object.kcur
– method order to be used on the next internal step.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetLastStep(void *ida_mem, realtype *hlast)
The function
IDAGetLastStep
returns the integration step size taken on the last internal step (if fromIDASolve()
), or the last value of the artificial step size \(h\) (if fromIDACalcIC()
). Arguments:
ida_mem
– pointer to the IDA solver object.hlast
– step size taken on the last internal step by IDA, or last artificial step size used inIDACalcIC()
, whichever was called last.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetCurrentStep(void *ida_mem, realtype *hcur)
The function
IDAGetCurrentStep
returns the integration step size to be attempted on the next internal step. Arguments:
ida_mem
– pointer to the IDA solver object.hcur
– step size to be attempted on the next internal step.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetActualInitStep(void *ida_mem, realtype *hinused)
The function
IDAGetActualInitStep
returns the value of the integration step size used on the first step. Arguments:
ida_mem
– pointer to the IDA solver object.hinused
– actual value of initial step size.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
Notes:
Even if the value of the initial integration step size was specified by the user through a call to
IDASetInitStep()
, this value might have been changed by IDA to ensure that the step size is within the prescribed bounds \((h_{min} \le h_0 \le h_{max})\), or to meet the local error test.

int IDAGetCurrentTime(void *ida_mem, realtype *tcur)
The function
IDAGetCurrentTime
returns the current internal time reached by the solver. Arguments:
ida_mem
– pointer to the IDA solver object.tcur
– current internal time reached.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetTolScaleFactor(void *ida_mem, realtype *tolsfac)
The function
IDAGetTolScaleFactor
returns a suggested factor by which the user’s tolerances should be scaled when too much accuracy has been requested for some internal step. Arguments:
ida_mem
– pointer to the IDA solver object.tolsfac
– suggested scaling factor for user tolerances.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetErrWeights(void *ida_mem, N_Vector eweight)
The function
IDAGetErrWeights
returns the solution error weights at the current time. These are the \(W_i\) given by (6.4) (or by the user’sIDAEwtFn
). Arguments:
ida_mem
– pointer to the IDA solver object.eweight
– solution error weights at the current time.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
Warning
The user must allocate space for
eweight
.

int IDAGetEstLocalErrors(void *ida_mem, N_Vector ele)
The function
IDAGetEstLocalErrors
returns the estimated local errors. Arguments:
ida_mem
– pointer to the IDA solver object.ele
– estimated local errors at the current time.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
Warning
The user must allocate space for
ele
. The values returned inele
are only valid ifIDASolve()
returned a nonnegative value.Note
The
ele
vector, togther with theeweight
vector fromIDAGetErrWeights()
, can be used to determine how the various components of the system contributed to the estimated local error test. Specifically, that error test uses the RMS norm of a vector whose components are the products of the components of these two vectors. Thus, for example, if there were recent error test failures, the components causing the failures are those with largest values for the products, denoted loosely aseweight[i]*ele[i]
.

int IDAGetIntegratorStats(void *ida_mem, long int *nsteps, long int *nrevals, long int *nlinsetups, long int *netfails, int *qlast, int *qcur, realtype *hinused, realtype *hlast, realtype *hcur, realtype *tcur)
The function
IDAGetIntegratorStats
returns the IDA integrator stats in one function call. Arguments:
ida_mem
– pointer to the IDA solver object.nsteps
– cumulative number of steps taken by IDA.nrevals
– cumulative number of calls to the user’sres
functions.nlinsetups
– cumulative number of calls made to the linear solver setup function.netfails
– cumulative number of error test failures.klast
– method order used on the last internal step.kcur
– method order to be used on the next internal step.hinused
– actual value of initial step size.hlast
– step sized taken on the last internal step.hcur
– step size to be attempted on the next internal step.tcur
– current internal time reached.
 Return value:
IDA_SUCCESS
– The optional output values have been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetNumNonlinSolvIters(void *ida_mem, long int *nniters)
The function
IDAGetNumNonlinSolvIters
returns the cumulative number of nonlinear iterations performed. Arguments:
ida_mem
– pointer to the IDA solver object.nniters
– number of nonlinear iterations performed.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_MEM_FAIL
– TheSUNNonlinearSolver
object isNULL
.

int IDAGetNumNonlinSolvConvFails(void *ida_mem, long int *nncfails)
The function
IDAGetNumNonlinSolvConvFails
returns the cumulative number of nonlinear convergence failures that have occurred. Arguments:
ida_mem
– pointer to the IDA solver object.nncfails
– number of nonlinear convergence failures.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetNonlinSolvStats(void *ida_mem, long int *nniters, long int *nncfails)
The function
IDAGetNonlinSolvStats
returns the IDA nonlinear solver statistics as a group. Arguments:
ida_mem
– pointer to the IDA solver object.nniters
– cumulative number of nonlinear iterations performed.nncfails
– cumulative number of nonlinear convergence failures.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_MEM_FAIL
– TheSUNNonlinearSolver
object isNULL
.

int IDAGetUserData(void *ida_mem, void **user_data)
The function
IDAGetUserData
returns the user data pointer provided toIDASetUserData()
. Arguments:
ida_mem
– pointer to the IDA memory block.user_data
– memory reference to a user data pointer.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
New in version 6.3.0.

int IDAPrintAllStats(void *ida_mem, FILE *outfile, SUNOutputFormat fmt)
The function
IDAPrintAllStats
outputs all of the integrator, nonlinear solver, linear solver, and other statistics. Arguments:
ida_mem
– pointer to the IDA memory block.outfile
– pointer to output file.fmt
– the output format:SUN_OUTPUTFORMAT_TABLE
– prints a table of valuesSUN_OUTPUTFORMAT_CSV
– prints a commaseparated list of key and value pairs e.g.,key1,value1,key2,value2,...
 Return value:
IDA_SUCCESS
– The output was successfully.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.IDA_ILL_INPUT
– An invalid formatting option was provided.
Note
The file
scripts/sundials_csv.py
provides python utility functions to read and output the data from a SUNDIALS CSV output file using the key and value pair format.New in version 6.2.0.

char *IDAGetReturnFlagName(long int flag)
The function
IDAGetReturnFlagName
returns the name of the IDA constant corresponding toflag
. Arguments:
flag
– the flag returned by a call to an IDA function
 Return value:
char*
– the flag name string
6.4.5.12.2. Initial condition calculation optional output functions

int IDAGetNumBacktrackOps(void *ida_mem, long int *nbacktr)
The function
IDAGetNumBacktrackOps
returns the number of backtrack operations done in the linesearch algorithm inIDACalcIC()
. Arguments:
ida_mem
– pointer to the IDA solver object.nbacktr
– the cumulative number of backtrack operations.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.

int IDAGetConsistentIC(void *ida_mem, N_Vector yy0_mod, N_Vector yp0_mod)
The function
IDAGetConsistentIC
returns the corrected initial conditions calculated byIDACalcIC()
. Arguments:
ida_mem
– pointer to the IDA solver object.yy0_mod
– consistent solution vector.yp0_mod
– consistent derivative vector.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_ILL_INPUT
– The function was not called before the first call toIDASolve()
.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
If the consistent solution vector or consistent derivative vector is not desired, pass
NULL
for the corresponding argument.
Warning
The user must allocate space for
yy0_mod
andyp0_mod
(if notNULL
).
6.4.5.12.3. Rootfinding optional output functions
There are two optional output functions associated with rootfinding.

int IDAGetRootInfo(void *ida_mem, int *rootsfound)
The function
IDAGetRootInfo
returns an array showing which functions were found to have a root. Arguments:
ida_mem
– pointer to the IDA solver object.rootsfound
– array of lengthnrtfn
with the indices of the user functions \(g_i\) found to have a root. For \(\mathtt{i} = 0, \ldots, \mathtt{nrtfn} 1\), \(\mathtt{rootsfound[i]} \ne 0\) if \(g_i\) has a root, and \(= 0\) if not.
 Return value:
IDA_SUCCESS
– The optional output values have been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
 Notes:
Note that, for the components \(g_i\) for which a root was found, the sign of
rootsfound[i]
indicates the direction of zerocrossing. A value of \(+1\) indicates that \(g_i\) is increasing, while a value of \(1\) indicates a decreasing \(g_i\).
Warning
The user must allocate memory for the vector
rootsfound
.

int IDAGetNumGEvals(void *ida_mem, long int *ngevals)
The function
IDAGetNumGEvals
returns the cumulative number of calls to the user root function \(g\). Arguments:
ida_mem
– pointer to the IDA solver object.ngevals
– number of calls to the user’s function \(g\) so far.
 Return value:
IDA_SUCCESS
– The optional output value has been successfully set.IDA_MEM_NULL
– Theida_mem
pointer isNULL
.
6.4.5.12.4. IDALS linear solver interface optional output functions
The following optional outputs are available from the IDALS modules:

int IDAGetLinWorkSpace(void *ida_mem, long int *lenrwLS, long int *leniwLS)
The function
IDAGetLinWorkSpace
returns the sizes of the real and integer workspaces used by the IDALS linear solver interface. Arguments:
ida_mem
– pointer to the IDA solver object.lenrwLS
– the number of real values in the IDALS workspace.leniwLS
– the number of integer values in the IDALS workspace.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
 Notes:
The workspace requirements reported by this routine correspond only to memory allocated within this interface and to memory allocated by the
SUNLinearSolver
object attached to it. The template Jacobian matrix allocated by the user outside of IDALS is not included in this report.
Warning
The previous routines
IDADlsGetWorkspace
andIDASpilsGetWorkspace
are now wrappers for this routine, and may still be used for backwardcompatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJacEvals(void *ida_mem, long int *njevals)
The function
IDAGetNumJacEvals
returns the cumulative number of calls to the IDALS Jacobian approximation function. Arguments:
ida_mem
– pointer to the IDA solver object.njevals
– the cumulative number of calls to the Jacobian function total so far.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDADlsGetNumJacEvals
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinResEvals(void *ida_mem, long int *nrevalsLS)
The function
IDAGetNumLinResEvals
returns the cumulative number of calls to the user residual function due to the finite difference Jacobian approximation or finite difference Jacobianvector product approximation. Arguments:
ida_mem
– pointer to the IDA solver object.nrevalsLS
– the cumulative number of calls to the user residual function.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
 Notes:
The value
nrevalsLS
is incremented only if one of the default internal difference quotient functions is used.
Warning
The previous routines
IDADlsGetNumRhsEvals
andIDASpilsGetNumRhsEvals
are now wrappers for this routine, and may still be used for backwardcompatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinIters(void *ida_mem, long int *nliters)
The function
IDAGetNumLinIters
returns the cumulative number of linear iterations. Arguments:
ida_mem
– pointer to the IDA solver object.nliters
– the current number of linear iterations.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumLinIters
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumLinConvFails(void *ida_mem, long int *nlcfails)
The function
IDAGetNumLinConvFails
returns the cumulative number of linear convergence failures. Arguments:
ida_mem
– pointer to the IDA solver object.nlcfails
– the current number of linear convergence failures.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumConvFails
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumPrecEvals(void *ida_mem, long int *npevals)
The function
IDAGetNumPrecEvals
returns the cumulative number of preconditioner evaluations, i.e., the number of calls made topsetup
. Arguments:
ida_mem
– pointer to the IDA solver object.npevals
– the cumulative number of calls topsetup
.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumPrecEvals
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumPrecSolves(void *ida_mem, long int *npsolves)
The function
IDAGetNumPrecSolves
returns the cumulative number of calls made to the preconditioner solve function,psolve
. Arguments:
ida_mem
– pointer to the IDA solver object.npsolves
– the cumulative number of calls topsolve
.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumPrecSolves
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJTSetupEvals(void *ida_mem, long int *njtsetup)
The function
IDAGetNumJTSetupEvals
returns the cumulative number of calls made to the Jacobianvector product setup functionjtsetup
. Arguments:
ida_mem
– pointer to the IDA solver object.njtsetup
– the current number of calls tojtsetup
.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumJTSetupEvals
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetNumJtimesEvals(void *ida_mem, long int *njvevals)
The function
IDAGetNumJtimesEvals
returns the cumulative number of calls made to the Jacobianvector product function,jtimes
. Arguments:
ida_mem
– pointer to the IDA solver object.njvevals
– the cumulative number of calls tojtimes
.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
Warning
The previous routine
IDASpilsGetNumJtimesEvals
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int IDAGetLastLinFlag(void *ida_mem, long int *lsflag)
The function
IDAGetLastLinFlag
returns the last return value from an IDALS routine. Arguments:
ida_mem
– pointer to the IDA solver object.lsflag
– the value of the last return flag from an IDALS function.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer isNULL
.IDALS_LMEM_NULL
– The IDALS linear solver has not been initialized.
 Notes:
If the IDALS setup function failed (i.e.,
IDASolve()
returnedIDA_LSETUP_FAIL
) when using the SUNLINSOL_DENSE or SUNLINSOL_BAND modules, then the value oflsflag
is equal to the column index (numbered from one) at which a zero diagonal element was encountered during the LU factorization of the (dense or banded) Jacobian matrix. If the IDALS setup function failed when using anotherSUNLinearSolver
object, thenlsflag
will beSUNLS_PSET_FAIL_UNREC
,SUNLS_ASET_FAIL_UNREC
, orSUNLS_PACKAGE_FAIL_UNREC
. If the IDALS solve function failed (IDASolve()
returnedIDA_LSOLVE_FAIL
),lsflag
contains the error return flag from theSUNLinearSolver
object, which will be one of:SUNLS_MEM_NULL
, indicating that theSUNLinearSolver
memory isNULL
;SUNLS_ATIMES_FAIL_UNREC
, indicating an unrecoverable failure in the \(J*v\) function;SUNLS_PSOLVE_FAIL_UNREC
, indicating that the preconditioner solve functionpsolve
failed unrecoverably;SUNLS_GS_FAIL
, indicating a failure in the GramSchmidt procedure (generated only in SPGMR or SPFGMR);SUNLS_QRSOL_FAIL
, indicating that the matrix \(R\) was found to be singular during the QR solve phase (SPGMR and SPFGMR only); orSUNLS_PACKAGE_FAIL_UNREC
, indicating an unrecoverable failure in an external iterative linear solver package.
Warning
The previous routines
IDADlsGetLastFlag
andIDASpilsGetLastFlag
are now wrappers for this routine, and may still be used for backwardcompatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

char *IDAGetLinReturnFlagName(long int lsflag)
The function
IDAGetLinReturnFlagName
returns the name of the IDALS constant corresponding tolsflag
. Arguments:
flag
– the flag returned by a call to an IDA function
 Return value:
char*
– the flag name string or if \(1 \leq \mathtt{lsflag} \leq N\) (LU factorization failed), this function returns “NONE”.
Warning
The previous routines
IDADlsGetReturnFlagName
andIDASpilsGetReturnFlagName
are now wrappers for this routine, and may still be used for backwardcompatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.
6.4.5.13. IDA reinitialization function
The function IDAReInit()
reinitializes the main IDA solver for the
solution of a new problem, where a prior call to IDAInit()
has been
made. The new problem must have the same size as the previous
one. IDAReInit()
performs the same input checking and initializations
that IDAInit()
does, but does no memory allocation, as it assumes that
the existing internal memory is sufficient for the new problem. A call to
IDAReInit()
deletes the solution history that was stored internally
during the previous integration. Following a successful call to
IDAReInit()
, call IDASolve()
again for the solution of the new
problem.
The use of IDAReInit()
requires that the maximum method order,
maxord
, is no larger for the new problem than for the problem specified in
the last call to IDAInit()
. In addition, the same N_Vector
module
set for the previous problem will be reused for the new problem.
If there are changes to the linear solver specifications, make the appropriate calls to either the linear solver objects themselves, or to the IDALS interface routines, as described in §6.4.5.5.
If there are changes to any optional inputs, make the appropriate IDASet***
calls, as described in §6.4.5.10.1. Otherwise,
all solver inputs set previously remain in effect.
One important use of the IDAReInit()
function is in the treating of jump
discontinuities in the residual function. Except in cases of fairly small jumps,
it is usually more efficient to stop at each point of discontinuity and restart
the integrator with a readjusted DAE model, using a call to IDAReInit()
.
To stop when the location of the discontinuity is known, simply make that
location a value of \(t_{\text{out}}\). To stop when the location of the
discontinuity is determined by the solution, use the rootfinding feature. In
either case, it is critical that the residual function not incorporate the
discontinuity, but rather have a smooth extention over the discontinuity, so
that the step across it (and subsequent rootfinding, if used) can be done
efficiently. Then use a switch within the residual function (communicated
through user_data
) that can be flipped between the stopping of the
integration and the restart, so that the restarted problem uses the new values
(which have jumped). Similar comments apply if there is to be a jump in the
dependent variable vector.

int IDAReInit(void *ida_mem, realtype t0, N_Vector y0, N_Vector yp0)
The function
IDAReInit
provides required problem specifications and reinitializes IDA. Arguments:
ida_mem
– pointer to the IDA solver object.t0
– is the initial value of \(t\).y0
– is the initial value of \(y\).yp0
– is the initial value of \(\dot{y}\).
 Return value:
IDA_SUCCESS
– The call to was successful.IDA_MEM_NULL
– The IDA solver object was not initialized through a previous call toIDACreate()
.IDA_NO_MALLOC
– Memory space for the IDA solver object was not allocated through a previous call toIDAInit()
.IDA_ILL_INPUT
– An input argument toIDAReInit()
has an illegal value.
 Notes:
If an error occurred,
IDAReInit()
also sends an error message to the error handler function.
6.4.6. Usersupplied functions
The usersupplied functions consist of one function defining the DAE residual, (optionally) a function that handles error and warning messages, (optionally) a function that provides the error weight vector, (optionally) one or two functions that provide Jacobianrelated information for the linear solver, and (optionally) one or two functions that define the preconditioner for use in any of the Krylov iteration algorithms.
6.4.6.1. DAE residual function
The user must provide a function of type IDAResFn
defined as follows:

typedef int (*IDAResFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, void *user_data)
This function computes the problem residual for given values of the independent variable \(t\), state vector \(y\), and derivative \(\dot{y}\).
 Arguments:
tt
– is the current value of the independent variable.yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the output residual vector \(F(t,y,\dot{y})\).user_data
– is a pointer to user data, the same as theuser_data
pointer parameter passed toIDASetUserData()
.
 Return value:
An
IDAResFn
function type should return a value of \(0\) if successful, a positive value if a recoverable error occurred (e.g.,yy
has an illegal value), or a negative value if a nonrecoverable error occurred. In the last case, the integrator halts. If a recoverable error occurred, the integrator will attempt to correct and retry. Notes:
A recoverable failure error return from the
IDAResFn
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, IDA will attempt to recover (possibly repeating the nonlinear solve, or reducing the step size) in order to avoid this recoverable error return.For efficiency reasons, the DAE residual function is not evaluated at the converged solution of the nonlinear solver. Therefore, in general, a recoverable error in that converged value cannot be corrected. (It may be detected when the residual function is called the first time during the following integration step, but a successful step cannot be undone.)
6.4.6.2. Error message handler function
As an alternative to the default behavior of directing error and warning
messages to the file pointed to by errfp
(see IDASetErrFile()
), the
user may provide a function of type IDAErrHandlerFn
to process any
such messages. The function type IDAErrHandlerFn
is defined as
follows:

typedef void (*IDAErrHandlerFn)(int error_code, const char *module, const char *function, char *msg, void *user_data)
This function processes error and warning messages from IDA and its submodules.
 Arguments:
error_code
– is the error code.module
– is the name of the IDA module reporting the error.function
– is the name of the function in which the error occurred.eH_data
– is a pointer to user data, the same as theeh_data
parameter passed toIDASetErrHandlerFn()
.
 Return value:
This function has no return value.
 Notes:
error_code
is negative for errors and positive (IDA_WARNING
) for warnings. If a function that returns a pointer to memory encounters an error, it setserror_code
to 0.
6.4.6.3. Error weight function

typedef int (*IDAEwtFn)(N_Vector y, N_Vector ewt, void *user_data)
This function computes the WRMS error weights for the vector \(y\).
 Arguments:
y
– is the value of the dependent variable vector at which the weight vector is to be computed.ewt
– is the output vector containing the error weights.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
0
– if it the error weights were successfully set.1
– if any error occured.
 Notes:
Allocation of memory for
ewt
is handled within IDA.
Warning
The error weight vector must have all components positive. It is the user’s responsiblity to perform this test and return 1 if it is not satisfied.
6.4.6.4. Rootfinding function
If a rootfinding problem is to be solved during the integration of the DAE
system, the user must supply a function of type IDARootFn
, defined
as follows:

typedef int (*IDARootFn)(realtype t, N_Vector y, N_Vector yp, realtype *gout, void *user_data)
This function computes a vectorvalued function \(g(t,y,\dot{y})\) such that the roots of the
nrtfn
components \(g_i(t,y,\dot{y})\) are to be found during the integration. Arguments:
t
– is the current value of the independent variable.y
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\), the \(t\text{derivative}\) of \(y\).gout
– is the output array, of lengthnrtfn
, with components \(g_i(t,y,\dot{y})\).user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
0
if successful or nonzero if an error occured (in which case the integration is halted andIDASolve()
retursIDA_RTFUNC_FAIL
). Notes:
Allocation of memory for
gout
is handled within IDA.
6.4.6.5. Jacobian construction (matrixbased linear solvers)
If a matrixbased linear solver module is used (i.e. a nonNULL
SUNMatrix
object was supplied to IDASetLinearSolver()
), the
user may provide a function of type IDALsJacFn
defined as follows:

typedef int (*IDALsJacFn)(realtype t, realtype c_j, N_Vector y, N_Vector yp, N_Vector r, SUNMatrix Jac, void *user_data, N_Vector tmp1, N_Vector tmp2, N_Vector tmp3)
This function computes the Jacobian matrix \(J\) of the DAE system (or an approximation to it), defined by (6.6).
 Arguments:
tt
– is the current value of the independent variable \(t\).cj
– is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the current value of the residual vector \(F(t,y,\dot{y})\).Jac
– is the output (approximate) Jacobian matrix (of typeSUNMatrix
), \(J = \dfrac{\partial{F}}{\partial{y}} + cj ~ \dfrac{\partial{F}}{\partial{\dot{y}}}\).user_data
 is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.tmp1
,tmp2
, andtmp3
– are pointers to memory allocated for variables of typeN_Vector
which can be used byIDALsJacFn
function as temporary storage or work space.
 Return value:
An
IDALsJacFn
should return \(0\) if successful, a positive value if a recoverable error occurred, or a negative value if a nonrecoverable error occurred.In the case of a recoverable eror return, the integrator will attempt to recover by reducing the stepsize, and hence changing \(\alpha\) in (6.6).
 Notes:
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 Chapter §10 for details).With direct linear solvers (i.e., linear solvers with type
SUNLINEARSOLVER_DIRECT
), the Jacobian matrix \(J(t,y,\dot{y})\) is zeroed out prior to calling the usersupplied Jacobian function so only nonzero elements need to be loaded intoJac
.With the default nonlinear solver (the native SUNDIALS Newton method), each call to the user’s
IDALsJacFn
function is preceded by a call to theIDAResFn
user function with the same (tt
,yy
,yp
) arguments. Thus the Jacobian function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual. In the case of a usersupplied or external nonlinear solver, this is also true if the residual function is evaluated prior to calling the linear solver setup function (see §12.1.4 for more information).If the user’s
IDALsJacFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These quantities may include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer toida_mem
touser_data
and then use theIDAGet*
functions described in §6.4.5.12.1. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.dense:
A usersupplied dense Jacobian function must load the
Neq
\(\times\)Neq
dense matrixJac
with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt
,yy
,yp
). The accessor macrosSM_ELEMENT_D
andSM_COLUMN_D
allow the user to read and write dense matrix elements without making explicit references to the underlying representation of theSUNMATRIX_DENSE
type.SM_ELEMENT_D(J, i, j)
references the (i
,j
)th element of the dense matrixJac
(withi
,j
\(= 0\ldots \texttt{N}1\)). This macro is meant for small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from \(1\) to \(N\), the Jacobian element \(J_{m,n}\) can be set using the statementSM_ELEMENT_D(J, m1, n1) =
\(J_{m,n}\). Alternatively,SM_COLUMN_D(J, j)
returns a pointer to the first element of thej
th column ofJac
(withj
\(= 0\ldots \texttt{N}1\)), and the elements of thej
th column can then be accessed using ordinary array indexing. Consequently, \(J_{m,n}\) can be loaded using the statementscol_n = SM_COLUMN_D(J, n1);
col_n[m1] =
\(J_{m,n}\). For large problems, it is more efficient to useSM_COLUMN_D
than to useSM_ELEMENT_D
. Note that both of these macros number rows and columns starting from \(0\). TheSUNMATRIX_DENSE
type and accessor macros are documented in §10.9.banded:
A usersupplied banded Jacobian function must load the
Neq
\(\times\)Neq
banded matrixJac
with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt
,yy
,yp
). The accessor macrosSM_ELEMENT_B
,SM_COLUMN_B
, andSM_COLUMN_ELEMENT_B
allow the user to read and write banded matrix elements without making specific references to the underlying representation of theSUNMATRIX_BAND
type.SM_ELEMENT_B(J, i, j)
references the (i
,j
)th element of the banded matrixJac
, counting from \(0\). This macro is meant for use in small problems for which efficiency of access is not a major concern. Thus, in terms of the indices \(m\) and \(n\) ranging from \(1\) to \(\texttt{N}\) with \((m,n)\) within the band defined bymupper
andmlower
, the Jacobian element \(J_{m,n}\) can be loaded using the statementSM_ELEMENT_B(J, m1, n1) =
\(J_{m,n}\). The elements within the band are those withmupper
\(\le\)mn
\(\le\)mlower
. Alternatively,SM_COLUMN_B(J, j)
returns a pointer to the diagonal element of thej
th column ofJac
, and if we assign this address torealtype *col_j
, then thei
th element of thej
th column is given bySM_COLUMN_ELEMENT_B(col_j, i, j)
, counting from \(0\). Thus, for \((m,n)\) within the band, \(J_{m,n}\) can be loaded by settingcol_n = SM_COLUMN_B(J, n1);
andSM_COLUMN_ELEMENT_B(col_n, m1, n1) =
\(J_{m,n}\). The elements of thej
th column can also be accessed via ordinary array indexing, but this approach requires knowledge of the underlying storage for a band matrix of typeSUNMATRIX_BAND
. The arraycol_n
can be indexed from \(\)mupper
tomlower
. For large problems, it is more efficient to useSM_COLUMN_B
andSM_COLUMN_ELEMENT_B
than to use theSM_ELEMENT_B
macro. As in the dense case, these macros all number rows and columns starting from \(0\). TheSUNMATRIX_BAND
type and accessor macros are documented in §10.12.sparse:
A usersupplied sparse Jacobian function must load the
Neq
\(\times\)Neq
compressedsparsecolumn or compressedsparserow matrixJac
with an approximation to the Jacobian matrix \(J(t,y,\dot{y})\) at the point (tt
,yy
,yp
). Storage forJac
already exists on entry to this function, although the user should ensure that sufficient space is allocated inJac
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. The amount of allocated space in aSUNMATRIX_SPARSE
object may be accessed using the macroSM_NNZ_S
or the routineSUNSparseMatrix_NNZ
. TheSUNMATRIX_SPARSE
type and accessor macros are documented in §10.14.
Warning
The previous function type
IDADlsJacFn
is identical toIDALsJacFn
, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.
6.4.6.6. Jacobianvector product (matrixfree linear solvers)
If a matrixfree linear solver is to be used (i.e., a NULL
valued
SUNMatrix
was supplied to IDASetLinearSolver()
), the user may
provide a function of type IDALsJacTimesVecFn
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 (*IDALsJacTimesVecFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector v, N_Vector Jv, realtype cj, void *user_data, N_Vector tmp1, N_Vector tmp2)
This function computes the product \(Jv\) of the DAE system Jacobian \(J\) (or an approximation to it) and a given vector
v
, where \(J\) is defined by (6.6). Arguments:
tt
– is the current value of the independent variable.yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the current value of the residual vector \(F(t,y,\dot{y})\).v
– is the vector by which the Jacobian must be multiplied to the right.Jv
– is the computed output vector.cj
– is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.tmp1
andtmp2
– are pointers to memory allocated for variables of typeN_Vector
which can be used byIDALsJacTimesVecFn
as temporary storage or work space.
 Return value:
The value returned by the Jacobiantimesvector function should be 0 if successful. A nonzero value indicates that a nonrecoverable error occurred.
 Notes:
This function must return a value of \(Jv\) that uses an approximation to the current value of \(J\), i.e. as evaluated at the current \((t,y,\dot{y})\).
If the user’s
IDALsJacTimesVecFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer toida_mem
touser_data
and then use theIDAGet*
functions described in §6.4.5.12.1. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.
Warning
The previous function type
IDASpilsJacTimesVecFn
is identical toIDALsJacTimesVecFn
, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.
6.4.6.7. Jacobianvector product setup (matrixfree linear solvers)
If the user’s Jacobianvector product function requires that any
Jacobianrelated data be preprocessed or evaluated, then this needs to be done
in a usersupplied function of type IDALsJacTimesSetupFn
, defined as
follows:

typedef int (*IDALsJacTimesSetupFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, ealtype cj, void *user_data);
This function setups any data needed by \(Jv\) product function (see
IDALsJacTimesVecFn
). Arguments:
tt
– is the current value of the independent variable.yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the current value of the residual vector \(F(t,y,\dot{y})\).cj
– is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
The value 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).
 Notes:
Each call to the Jacobianvector product setup function is preceded by a call to the
IDAResFn
user function with the same \((t,y,\dot{y})\) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual.If the user’s
IDALsJacTimesVecFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer toida_mem
touser_data
and then use theIDAGet*
functions described in §6.4.5.12.1. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.
Warning
The previous function type
IDASpilsJacTimesSetupFn
is identical toIDALsJacTimesSetupFn
, and may still be used for backwardcompatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.
6.4.6.8. Preconditioner solve (iterative linear solvers)
If a usersupplied preconditioner is to be used with a SUNLinearSolver
solver module, then the user must provide a function to solve the linear system
\(Pz = r\) where \(P\) is a left preconditioner matrix which
approximates (at least crudely) the Jacobian matrix \(J =
\dfrac{\partial{F}}{\partial{y}} + cj ~ \dfrac{\partial{F}}{\partial{\dot{y}}}\). This function
must be of type IDALsPrecSolveFn
, defined as follows:

typedef int (*IDALsPrecSolveFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, N_Vector rvec, N_Vector zvec, realtype cj, realtype delta, void *user_data)
This function solves the preconditioning system \(Pz = r\).
 Arguments:
tt
– is the current value of the independent variable.yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the current value of the residual vector \(F(t,y,\dot{y})\).rvec
– is the righthand side vector \(r\) of the linear system to be solved.zvec
– is the computed output vector.cj
– is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).delta
– is an input tolerance to be used if an iterative method is employed in the solution. In that case, the residual vector \(Res = r  P z\) of the system should be made less thandelta
in weighted \(l_2\) norm, i.e., \(\sqrt{\displaystyle\sum_i (Res_i \cdot ewt_i)^2 } <\)delta
. To obtain theN_Vector
ewt
, callIDAGetErrWeights()
.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
The value returned by the preconditioner solve 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).
6.4.6.9. Preconditioner setup (iterative linear solvers)
If the user’s preconditioner requires that any Jacobianrelated data be
evaluated or preprocessed, then this needs to be done in a usersupplied
function of type IDALsPrecSetupFn
, defined as follows:

typedef int (*IDALsPrecSetupFn)(realtype tt, N_Vector yy, N_Vector yp, N_Vector rr, realtype cj, void *user_data)
This function solves the preconditioning system \(Pz = r\).
 Arguments:
tt
– is the current value of the independent variable.yy
– is the current value of the dependent variable vector, \(y(t)\).yp
– is the current value of \(\dot{y}(t)\).rr
– is the current value of the residual vector \(F(t,y,\dot{y})\).cj
– is the scalar in the system Jacobian, proportional to the inverse of the step size (\(\alpha\) in (6.6)).user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
The value returned by the preconditioner 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).
 Notes:
With the default nonlinear solver (the native SUNDIALS Newton method), each call to the preconditioner setup function is preceded by a call to the
IDAResFn
user function with the same \((t,y,\dot{y})\) arguments. Thus the preconditioner setup function can use any auxiliary data that is computed and saved during the evaluation of the DAE residual. In the case of a usersupplied or external nonlinear solver, this is also true if the residual function is evaluated prior to calling the linear solver setup function (see §12.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 nonlinear solver.
If the user’s
IDALsPrecSetupFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer toida_mem
touser_data
and then use theIDAGet*
functions described in §6.4.5.12.1. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.
6.4.7. Preconditioner modules
A principal reason for using a parallel DAE solver such as IDA lies in the solution of partial differential equations (PDEs). Moreover, the use of a Krylov iterative method for the solution of many such problems is motivated by the nature of the underlying linear system of equations (6.5) that must be solved at each time step. The linear algebraic system is large, sparse, and structured. However, if a Krylov iterative method is to be effective in this setting, then a nontrivial preconditioner needs to be used. Otherwise, the rate of convergence of the Krylov iterative method is usually unacceptably slow. Unfortunately, an effective preconditioner tends to be problemspecific.
However, we have developed one type of preconditioner that treats a rather broad class of PDEbased problems. It has been successfully used for several realistic, largescale problems [65] and is included in a software module within the IDA package. This module works with the parallel vector module NVECTOR_PARALLEL and generates a preconditioner that is a blockdiagonal matrix with each block being a band matrix. The blocks need not have the same number of super and subdiagonals, and these numbers may vary from block to block. This BandBlockDiagonal Preconditioner module is called IDABBDPRE.
One way to envision these preconditioners is to think of the domain of the computational PDE problem as being subdivided into \(M\) nonoverlapping subdomains. Each of these subdomains is then assigned to one of the \(M\) processors to be used to solve the DAE system. The basic idea is to isolate the preconditioning so that it is local to each processor, and also to use a (possibly cheaper) approximate residual function. This requires the definition of a new function \(G(t,y,\dot{y})\) which approximates the function \(F(t, y, \dot{y})\) in the definition of the DAE system (6.1). However, the user may set \(G = F\). Corresponding to the domain decomposition, there is a decomposition of the solution vectors \(y\) and \(\dot{y}\) into \(M\) disjoint blocks \(y_m\) and \(\dot{y}_m\), and a decomposition of \(G\) into blocks \(G_m\). The block \(G_m\) depends on \(y_m\) and \(\dot{y}_m\), and also on components of \(y_{m'}\) and \(\dot{y}_{m'}\) associated with neighboring subdomains (socalled ghostcell data). Let \(\bar{y}_m\) and \(\bar{\dot{y}}_m\) denote \(y_m\) and \(\dot{y}_m\) (respectively) augmented with those other components on which \(G_m\) depends. Then we have
and each of the blocks \(G_m(t,\bar{y}_m,\bar{\dot{y}}_m)\) is uncoupled from the others.
The preconditioner associated with this decomposition has the form
where
This matrix is taken to be banded, with upper and lower halfbandwidths mudq
and mldq
defined as the number of nonzero diagonals above and below the
main diagonal, respectively. The difference quotient approximation is computed
using mudq
\(+\) mldq
\(+ 2\) evaluations of \(G_m\), but
only a matrix of bandwidth mukeep
\(+\) mlkeep
\(+ 1\) is
retained.
Neither pair of parameters need be the true halfbandwidths of the Jacobians of
the local block of \(G\), if smaller values provide a more efficient
preconditioner. Such an efficiency gain may occur if the couplings in the DAE
system outside a certain bandwidth are considerably weaker than those within the
band. Reducing mukeep
and mlkeep
while keeping mudq
and mldq
at
their true values, discards the elements outside the narrower band. Reducing
both pairs has the additional effect of lumping the outer Jacobian elements into
the computed elements within the band, and requires more caution and
experimentation.
The solution of the complete linear system
reduces to solving each of the equations
and this is done by banded LU factorization of \(P_m\) followed by a banded backsolve.
Similar blockdiagonal preconditioners could be considered with different treatment of the blocks \(P_m\). For example, incomplete LU factorization or an iterative method could be used instead of banded LU factorization.
6.4.7.1. A parallel bandblockdiagonal preconditioner module
The IDABBDPRE module calls two userprovided functions to construct \(P\): a
required function Gres
(of type IDABBDLocalFn
) which approximates
the residual function \(G(t,y,\dot{y}) \approx F(t,y,\dot{y})\) and which is
computed locally, and an optional function Gcomm
(of type
IDABBDCommFn
) which performs all interprocess communication necessary
to evaluate the approximate residual \(G\). These are in addition to the
usersupplied residual function res
. Both functions take as input the same
pointer user_data
as passed by the user to IDASetUserData()
and
passed to the user’s function res
. The user is responsible for providing
space (presumably within user_data
) for components of yy
and yp
that
are communicated by Gcomm
from the other processors, and that are then used
by Gres
, which should not do any communication.

typedef int (*IDABBDLocalFn)(sunindextype Nlocal, realtype tt, N_Vector yy, N_Vector yp, N_Vector gval, void *user_data)
This
Gres
function computes \(G(t,y,\dot{y})\). It loads the vectorgval
as a function oftt
,yy
, andyp
. Arguments:
Nlocal
– is the local vector length.tt
– is the value of the independent variable.yy
– is the dependent variable.yp
– is the derivative of the dependent variable.gval
– is the output vector.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
An
IDABBDLocalFn
function type should return 0 to indicate success, 1 for a recoverable error, or 1 for a nonrecoverable error. Notes:
This function must assume that all interprocessor communication of data needed to calculate
gval
has already been done, and this data is accessible withinuser_data
.The case where \(G\) is mathematically identical to \(F\) is allowed.

typedef int (*IDABBDCommFn)(sunindextype Nlocal, realtype tt, N_Vector yy, N_Vector yp, void *user_data)
This
Gcomm
function performs all interprocessor communications necessary for the execution of theGres
function above, using the input vectorsyy
andyp
. Arguments:
Nlocal
– is the local vector length.tt
– is the value of the independent variable.yy
– is the dependent variable.yp
– is the derivative of the dependent variable.gval
– is the output vector.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toIDASetUserData()
.
 Return value:
An
IDABBDCommFn
function type should return 0 to indicate success, 1 for a recoverable error, or 1 for a nonrecoverable error. Notes:
The
Gcomm
function is expected to save communicated data in space defined within the structureuser_data
.Each call to the
Gcomm
function is preceded by a call to the residual functionres
with the same \((t,y,\dot{y})\) arguments. ThusGcomm
can omit any communications done byres
if relevant to the evaluation ofGres
. If all necessary communication was done inres
, thenGcomm = NULL
can be passed in the call toIDABBDPrecInit()
.
Besides the header files required for the integration of the DAE problem (see
§6.4.3), to use the IDABBDPRE module, the main program
must include the header file ida_bbdpre.h
which declares the needed function
prototypes.
The following is a summary of the usage of this module and describes the sequence of calls in the user main program. Steps that are unchanged from the user main program presented in §6.4.4 are not bold.
Initialize parallel or multithreaded environment (if appropriate)
Create the vector of initial values
Create matrix object (if appropriate)
Create linear solver object (if appropriate)
When creating the iterative linear solver object, specify the use of left preconditioning (
SUN_PREC_LEFT
) as IDA only supports left preconditioning.Create nonlinear solver object (if appropriate)
Create IDA object
Initialize IDA solver
Specify integration tolerances
Attach the linear solver (if appropriate)
Set linear solver optional inputs (if appropriate)
Note that the user should not overwrite the preconditioner setup function or solve function through calls to
IDASetPreconditioner()
optional input function.Initialize the IDABBDPRE preconditioner module
Call
IDABBDPrecInit()
to allocate memory and initialize the internal preconditioner data. The last two arguments ofIDABBDPrecInit()
are the two usersupplied functions described above.Attach nonlinear solver module (if appropriate)
Set nonlinear solver optional inputs (if appropriate)
Specify rootfinding problem (optional)
Set optional inputs
Advance solution in time
Get optional outputs
Additional optional outputs associated with IDABBDPRE are available by way of two routines described below,
IDABBDPrecGetWorkSpace()
andIDABBDPrecGetNumGfnEvals()
.Deallocate memory
Finalize MPI, if used
The usercallable functions that initialize or reinitialize the IDABBDPRE preconditioner module are described next.

int IDABBDPrecInit(void *ida_mem, sunindextype Nlocal, sunindextype mudq, sunindextype mldq, sunindextype mukeep, sunindextype mlkeep, realtype dq_rel_yy, IDABBDLocalFn Gres, IDABBDCommFn Gcomm);
The function
IDABBDPrecInit
initializes and allocates (internal) memory for the IDABBDPRE preconditioner. Arguments:
ida_mem
– pointer to the IDA solver object.Nlocal
– local vector dimension.mudq
– upper halfbandwidth to be used in the differencequotient Jacobian approximation.mldq
– lower halfbandwidth to be used in the differencequotient Jacobian approximation.mukeep
– upper halfbandwidth of the retained banded approximate Jacobian block.mlkeep
– lower halfbandwidth of the retained banded approximate Jacobian block.dq_rel_yy
– the relative increment in components ofy
used in the difference quotient approximations. The default is \(\mathtt{dq\_rel\_yy} = \sqrt{\text{unit roundoff}}\) , which can be specified by passing \(\mathtt{dq\_rel\_yy} = 0.0\).Gres
– the function which computes the local residual approximation \(G(t,y,\dot{y})\).Gcomm
– the optional function which performs all interprocess communication required for the computation of \(G(t,y,\dot{y})\).
 Return value:
IDALS_SUCCESS
– The call was successful.IDALS_MEM_NULL
– Theida_mem
pointer wasNULL
.IDALS_MEM_FAIL
– A memory allocation request has failed.IDALS_LMEM_NULL
– An IDALS linear solver memory was not attached.IDALS_ILL_INPUT
– The supplied vector implementation was not compatible with the block band preconditioner.
 Notes:
If one of the halfbandwidths
mudq
ormldq
to be used in the differencequotient calculation of the approximate Jacobian is negative or exceeds the valueNlocal1
, it is replaced by 0 orNlocal1
accordingly. The halfbandwidthsmudq
andmldq
need not be the true halfbandwidths of the Jacobian of the local block of \(G\), when smaller values may provide a greater efficiency. Also, the halfbandwidthsmukeep
andmlkeep
of the retained banded approximate Jacobian block may be even smaller, to reduce storage and computation costs further. For all four halfbandwidths, the values need not be the same on every processor.
The IDABBDPRE module also provides a reinitialization function to allow for a
sequence of problems of the same size, with the same linear solver choice,
provided there is no change in local_N
, mukeep
, or mlkeep
. After
solving one problem, and after calling IDAReInit()
to reinitialize IDA
for a subsequent problem, a call to IDABBDPrecReInit()
can be made to
change any of the following: the halfbandwidths mudq
and mldq
used in
the differencequotient Jacobian approximations, the relative increment
dq_rel_yy
, or one of the usersupplied functions Gres
and Gcomm
. If
there is a change in any of the linear solver inputs, an additional call to the
“Set”routines provided by the SUNLinearSolver
object, and/or one or more of
the corresponding IDASet***
functions, must also be made (in the proper
order).

int IDABBDPrecReInit(void *ida_mem, sunindextype mudq, sunindextype mldq, realtype dq_rel_yy)
The function
IDABBDPrecReInit
reinitializes the IDABBDPRE preconditioner. Arguments:
ida_mem
– pointer to the IDA solver object.mudq
– upper halfbandwidth to be used in the differencequotient Jacobian approximation.Mldq
– lower halfbandwidth to be used in the differencequotient Jacobian approximation.dq_rel_yy
– the relative increment in components ofy
used in the difference quotient approximations. The default is \(\mathtt{dq\_rel\_yy} = \sqrt{\text{unit roundoff}}\) , which can be specified by passing \(\mathtt{dq\_rel\_yy} = 0.0\).
 Return value:
IDALS_SUCCESS
– The call was successful.IDALS_MEM_NULL
– Theida_mem
pointer wasNULL
.IDALS_LMEM_NULL
– An IDALS linear solver memory was not attached.IDALS_PMEM_NULL
– The functionIDABBDPrecInit()
was not previously called.
 Notes:
If one of the halfbandwidths
mudq
ormldq
is negative or exceeds the valueNlocal  1
, it is replaced by 0 orNlocal  1
, accordingly.
The following two optional output functions are available for use with the IDABBDPRE module:

int IDABBDPrecGetWorkSpace(void *ida_mem, long int *lenrwBBDP, long int *leniwBBDP)
The function
IDABBDPrecGetWorkSpace
returns the local sizes of the IDABBDPRE real and integer workspaces. Arguments:
ida_mem
– pointer to the IDA solver object.lenrwBBDP
– local number of real values in the IDABBDPRE workspace.leniwBBDP
– local number of integer values in the IDABBDPRE workspace.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer wasNULL
.IDALS_PMEM_NULL
– The IDABBDPRE preconditioner has not been initialized.
 Notes:
The workspace requirements reported by this routine correspond only to memory allocated within the IDABBDPRE module (the banded matrix approximation, banded
SUNLinearSolver
object, temporary vectors). These values are local to each process. The workspaces referred to here exist in addition to those given by the corresponding functionIDAGetLinWorkSpace()
.

int IDABBDPrecGetNumGfnEvals(void *ida_mem, long int *ngevalsBBDP)
The function
IDABBDPrecGetNumGfnEvals
returns the cumulative number of calls to the userGres
function due to the finite difference approximation of the Jacobian blocks used within IDABBDPRE’s preconditioner setup function. Arguments:
ida_mem
– pointer to the IDA solver object.ngevalsBBDP
– the cumulative number of calls to the userGres
function.
 Return value:
IDALS_SUCCESS
– The optional output value has been successfully set.IDALS_MEM_NULL
– Theida_mem
pointer wasNULL
.IDALS_PMEM_NULL
– The IDABBDPRE preconditioner has not been initialized.
In addition to the ngevalsBBDP
evaluations of Gres
, the costs associated
with IDABBDPRE also includes nlinsetups
LU factorizations, nlinsetups
calls to Gcomm
, npsolves
banded backsolve calls, and nrevalsLS
residual function evaluations, where nlinsetups
is an optional IDA output
(see §6.4.5.12.1), and npsolves
and
nrevalsLS
are linear solver optional outputs (see
§6.4.5.12.4).