8.4. Using KINSOL for the Solution of Nonlinear Systems
This section is concerned with the use of KINSOL for the solution of nonlinear systems.
The following sections treat the header files and the layout of the user’s main program, and provide descriptions of the KINSOL usercallable functions and usersupplied functions. The sample programs described in the companion document [38] 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 KINSOL package.
KINSOL uses various constants for both input and output. These are defined as needed in this chapter, but for convenience are also listed separately in §8.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 KINBBDPRE
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.
8.4.1. Access to library and header files
At this point, it is assumed that the installation of KINSOL, 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 KINSOL. The relevant library files are
<libdir>/libsundials_kinsol.<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>/kinsol
<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).
8.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
8.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).
8.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).
8.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
).
8.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.
8.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:
kinsol/kinsol.h
the main header file for kinsol, which defines the types and various constants, and includes function prototypes. This includes the header file for KINLS,kinsol/kinsol_ls.h
.
Note that kinsol.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 §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 Newton or Picard nonlinear solver that requires the solution of a
linear system, then a linear solver module header file will be required. 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 kinFoodWeb_kry_p
(see [38]),
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.
8.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 solution of a nonlinear system problem.. Most of the steps are independent
of the N_Vector
, SUNMatrix
, and SUNLinearSolver
implementations
used. For the steps that are not, refer to §9,
§10, and §11 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.Set the problem dimensions etc.
This generally includes the problem size
N
, and may include the local vector length Nlocal.Create the vector with the initial guess
Construct an
N_Vector
of initial guess 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.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 KINSOL object
Call
KINCreate()
to create the KINSOL solver object.Initialize KINSOL solver
Call
KINInit()
to allocate internal memory.Attach the linear solver (if appropriate)
If a linear solver was created above, initialize the KINLS linear solver interface by attaching the linear solver object (and matrix object, if applicable) with
KINSetLinearSolver()
.Set linear solver optional inputs (if appropriate)
See Table 8.1 for KINLS optional inputs and Chapter §11 for linear solver specific optional inputs.
Set optional inputs
Call
KINSet***
functions to change any optional inputs that control the behavior of KINSOL from their default values. See §8.4.5.4 for details.Solve problem
Call
ier = KINSol(...)
to solve the nonlinear problem for a given initial guess.See
KINSol()
for details.Get optional outputs
Call
KINGet***
functions to obtain optional output. See §8.4.5.5 for details.Deallocate memory
Upon completion of the integration call the following, as necessary, to free any objects or memory allocated 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
KINFree()
to free the memory allocated by KINSOL.Call
SUNContext_Free()
to free theSUNContext
object
Finalize MPI, if used
Call
MPI_Finalize
to terminate MPI.
8.4.5. Usercallable functions
This section describes the KINSOL functions that are called by the user to setup and then solve an IVP. Some of these are required. However, starting with §8.4.5.4, the functions listed involve optional inputs/outputs or restarting, and those paragraphs may be skipped for a casual use of KINSOL. In any case, refer to §8.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 §8.4.5.4).
8.4.5.1. KINSOL initialization and deallocation functions

void KINCreate(SUNContext sunctx)
The function
KINCreate()
instantiates a KINSOL solver object. Arguments:
sunctx
– theSUNContext
object (see §2.1)
 Return value:
void

int KINInit(void *kin_mem, KINSysFn func, N_Vector tmpl)
The function
KINInit()
specifies the problemdefining function, allocates internal memory, and initializes KINSOL. Arguments:
kin_mem
– pointer to the KINSOL memory block returned byKINCreate()
.func
– is the CC function which computes the system function \(F(u)\) (or \(G(u)\) for fixedpoint iteration) in the nonlinear problem. This function has the formfunc(u, fval, user_data)
. (For full details see §8.4.6.1).tmpl
– is anyN_Vector
(e.g. the initial guess vectoru
) which is used as a template to create (by cloning) necessary vectors inkin_mem
.
 Return value:
KIN_SUCCESS
– The call toKINInit()
was successful.KIN_MEM_NULL
– The KINSOL memory block was not initialized through a previous call toKINCreate()
.KIN_MEM_FAIL
– A memory allocation request has failed.KIN_ILL_INPUT
– An input argument toKINInit()
has an illegal value.
 Notes:
If an error occurred,
KINInit()
sends an error message to the error handler function.

void KINFree(void **kin_mem)
The function
KINFree()
frees the pointer allocated by a previous call toKINCreate()
. Arguments:
kin_mem
– pointer to the KINSOL solver object.
 Return value:
void
8.4.5.2. Linear solver specification functions
As previously explained, Newton and Picard iterations require the solution of
linear systems of the form \(J\delta = F\). Solution of these linear
systems is handled using the KINLS linear solver interface. This interface
supports all valid SUNLinearSolver
modules. Here, matrixbased
SUNLinearSolver
modules utilize SUNMatrix
objects to store the Jacobian
matrix \(J = F'(u)\) and factorizations used throughout
the solution process. Conversely, matrixfree SUNLinearSolver
modules
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. However, only right preconditioning is supported within KINLS. If preconditioning is done, usersupplied functions define the linear operator corresponding to a right preconditioner matrix \(P\), which should approximate the system Jacobian matrix \(J\). For the specification of a preconditioner, see the iterative linear solver sections in §8.4.5.4 and §8.4.6. A preconditioner matrix \(P\) must approximate the Jacobian \(J\), at least crudely.
To specify a generic linear solver to KINSOL, after the call to KINCreate()
but before any calls to KINSol()
, the user’s program must create the
appropriate SUNLinearSolver
object and call the function
KINSetLinearSolver()
, as documented below. To create the SUNLinearSolver
object, the user may call one of the SUNDIALSpackaged SUNLinearSolver
module constructor routines via a call of the form
SUNLinearSolver LS = SUNLinSol_*(...);
For a current list of such constructor routines see §11.
Alternately, a usersupplied SUNLinearSolver
module 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
module in question, as
described in Chapters §10 and §11.
Once this solver object has been constructed, the user should attach it to
KINSOL via a call to KINSetLinearSolver()
. The first argument passed to this
function is the KINSOL memory pointer returned by KINCreate()
; the second
argument is the desired SUNLinearSolver
object to use for solving Newton or
Picard 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 KINLS linear solver interface, linking it to the main
KINSOL solver, and allows the user to specify additional parameters and routines
pertinent to their choice of linear solver.

int KINSetLinearSolver(void *kin_mem, SUNLinearSolver LS, SUNMatrix J)
The function
KINSetLinearSolver()
attaches a genericSUNLinSol
objectLS
and corresponding template JacobianSUNMatrix
objectJ
(if applicable) to KINSOL, initializing the KINLS linear solver interface. Arguments:
kin_mem
– pointer to the KINSOL memory block.LS
– SUNLINSOL object to use for solving Newton linear systems.J
– SUNMATRIX object for used as a template for the Jacobian (orNULL
if not applicable).
 Return value:
KINLS_SUCCESS
– The KINLS initialization was successful.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_ILL_INPUT
– The KINLS interface is not compatible with theLS
orJ
input objects or is incompatible with the current NVECTOR module.KINLS_SUNLS_FAIL
– A call to theLS
object failed.KINLS_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
type in Chapter §10 for further information).The previous routines
KINDlsSetLinearSolver()
andKINSpilsSetLinearSolver()
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.
8.4.5.3. KINSOL solver function
This is the central step in the solution process, the call to solve the nonlinear algebraic system.

int KINSol(void *kin_mem, N_Vector u, int strategy, N_Vector u_scale, N_Vector f_scale)
The function
KINSol()
computes an approximate solution to the nonlinear system. Arguments:
kin_mem
– pointer to the KINSOL memory block.u
– vector set to initial guess by user before callingKINSol()
, but which upon return contains an approximate solution of the nonlinear system \(F(u) = 0\).strategy
– strategy used to solve the nonlinear system. It must be of the following:KIN_NONE
basic Newton iterationKIN_LINESEARCH
Newton with globalizationKIN_FP
fixedpoint iteration with Anderson Acceleration (no linear solver needed)KIN_PICARD
Picard iteration with Anderson Acceleration (uses a linear solver)
u_scale
– vector containing diagonal elements of scaling matrix \(D_u\) for vectoru
chosen so that the components of \(D_u\ u\) (as a matrix multiplication) all have roughly the same magnitude whenu
is close to a root of \(F(u)\).f_scale
– vector containing diagonal elements of scaling matrix \(D_F\) for \(F(u)\) chosen so that the components of \(D_F\ F(u)\) (as a matrix multiplication) all have roughly the same magnitude whenu
is not too near a root of \(F(u)\). In the case of a fixedpoint iteration, consider \(F(u) = G(u)  u\).
 Return value:
KIN_SUCCESS
–KINSol()
succeeded; the scaled norm of \(F(u)\) is less thanfnormtol
.KIN_INITIAL_GUESS_OK
– The guessu
\(=u_0\) satisfied the system \(F(u)=0\) within the tolerances specified (the scaled norm of \(F(u_0)\) is less than0.01*fnormtol
).KIN_STEP_LT_STPTOL
– KINSOL stopped based on scaled step length. This means that the current iterate may be an approximate solution of the given nonlinear system, but it is also quite possible that the algorithm is “stalled” (making insufficient progress) near an invalid solution, or that the scalarscsteptol
is too large (seeKINSetScaledStepTol()
in §8.4.5.4 to changescsteptol
from its default value).KIN_MEM_NULL
– The KINSOL memory block pointer wasNULL
.KIN_ILL_INPUT
– An input parameter was invalid.KIN_NO_MALLOC
– The KINSOL memory was not allocated by a call toKINCreate()
.KIN_MEM_FAIL
– A memory allocation failed.KIN_LINESEARCH_NONCONV
– The line search algorithm was unable to find an iterate sufficiently distinct from the current iterate, or could not find an iterate satisfying the sufficient decrease condition. Failure to satisfy the sufficient decrease condition could mean the current iterate is “close” to an approximate solution of the given nonlinear system, the difference approximation of the matrixvector product \(J(u)\ v\) is inaccurate, or the real scalarscsteptol
is too large.KIN_MAXITER_REACHED
– The maximum number of nonlinear iterations has been reached.KIN_MXNEWT_5X_EXCEEDED
– Five consecutive steps have been taken that satisfy the inequality \(\D_u p\_{L2} > 0.99\ \texttt{mxnewtstep}\) , where \(p\) denotes the current step andmxnewtstep
is a scalar upper bound on the scaled step length. Such a failure may mean that \(\D_F F(u)\_{L2}\) asymptotes from above to a positive value, or the real scalarmxnewtstep
is too small.KIN_LINESEARCH_BCFAIL
– The line search algorithm was unable to satisfy the “betacondition” forMXNBCF+1
nonlinear iterations (not necessarily consecutive), which may indicate the algorithm is making poor progress.KIN_LINSOLV_NO_RECOVERY
– The usersupplied routinepsolve
encountered a recoverable error, but the preconditioner is already current.KIN_LINIT_FAIL
– The KINLS initialization routine (linit
) encountered an error.KIN_LSETUP_FAIL
– The KINLS setup routine (lsetup
) encountered an error; e.g., the usersupplied routinepset
(used to set up the preconditioner data) encountered an unrecoverable error.KIN_LSOLVE_FAIL
– The KINLS solve routine (lsolve
) encountered an error; e.g., the usersupplied routinepsolve
(used to to solve the preconditioned linear system) encountered an unrecoverable error.KIN_SYSFUNC_FAIL
– The system function failed in an unrecoverable manner.KIN_FIRST_SYSFUNC_ERR
– The system function failed recoverably at the first call.KIN_REPTD_SYSFUNC_ERR
– The system function had repeated recoverable errors. No recovery is possible.
 Notes:
The components of vectors
u_scale
andf_scale
should be strictly positive.KIN_SUCCESS=0
,KIN_INITIAL_GUESS_OK=1
, andKIN_STEP_LT_STPTOL=2
. All remaining return values are negative and therefore a testflag
\(< 0\) will trap allKINSol()
failures.
8.4.5.4. Optional input functions
There are numerous optional input parameters that control the behavior of the KINSOL solver. KINSOL provides functions that can be used to change these from their default values. Table 8.1 lists all optional input functions in KINSOL which are then described in detail in the remainder of this section, beginning with those for the main KINSOL solver and continuing with those for the KINLS linear solver interface.
We note that, on error return, all of these functions also send an error message
to the error handler function. We also note that all error return values are
negative, so a test retval
\(<0\) will catch any error.
Optional input 
Function name 
Default 

KINSOL main solver 

Error handler function 
internal fn. 

Pointer to an error file 


Info handler function 
internal fn. 

Pointer to an info file 


Data for problemdefining function 


Verbosity level of output 
0 

Max. number of nonlinear iterations 
200 

No initial matrix setup 


No residual monitoring 


Max. iterations without matrix setup 
10 

Max. iterations without residual check 
5 

Form of \(\eta\) coefficient 


Constant value of \(\eta\) 
0.1 

Values of \(\gamma\) and \(\alpha\) 
0.9 and 2.0 

Values of \(\omega_{min}\) and \(\omega_{max}\) 
0.00001 and 0.9 

Constant value of \(\omega\) 
0.9 

Lower bound on \(\epsilon\) 


Max. scaled length of Newton step 
\(1000D_u u_0_2\) 

Max. number of \(\beta\)condition failures 
10 

Rel. error for D.Q. \(Jv\) 
\(\sqrt{\text{uround}}\) 

Functionnorm stopping tolerance 
uround\(^{1/3}\) 

Scaledstep stopping tolerance 
\(\text{uround}^{2/3}\) 

Inequality constraints on solution 


Nonlinear system function 
none 

Return the newest fixed point iteration 


Fixed point/Picard damping parameter 
1.0 

Anderson Acceleration subspace size 
0 

Anderson Acceleration damping parameter 
1.0 

Anderson Acceleration delay 
0 

Anderson Acceleration orthogonalization routine 


KINLS linear solver interface 

Jacobian function 
DQ 

Preconditioner functions and data 


Jacobiantimesvector function and data 
internal DQ, 

Jacobiantimesvector system function 


int KINSetErrFile(void *kin_mem, FILE *errfp)
The function
KINSetErrFile()
specifies the pointer to the file where all KINSOL messages should be directed when the default KINSOL error handler function is used. Arguments:
kin_mem
– pointer to the KINSOL memory block.errfp
– pointer to output file.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default value for
errfp
isstderr
.Passing a value of
NULL
disables all future error message output (except for the case in which the KINSOL memory pointer isNULL
). This use ofKINSetErrFile()
is strongly discouraged.
Warning
If
KINSetErrFile()
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 KINSetErrHandlerFn(void *kin_mem, KINErrHandlerFn ehfun, void *eh_data)
The function
KINSetErrHandlerFn()
specifies the optional userdefined function to be used in handling error messages. Arguments:
kin_mem
– pointer to the KINSOL memory block.ehfun
– is the user’s CC error handler function (see §8.4.6.2).eh_data
– pointer to user data passed toehfun
every time it is called.
 Return value:
KIN_SUCCESS
– The functionehfun
and data pointereh_data
have been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default internal error handler function directs error messages to the file specified by the file pointer
errfp
(seeKINSetErrFile()
above).Error messages indicating that the KINSOL solver memory is
NULL
will always be directed tostderr
.

int KINSetInfoFile(void *kin_mem, FILE *infofp)
The function
KINSetInfoFile()
specifies the pointer to the file where all informative (nonerror) messages should be directed. Arguments:
kin_mem
– pointer to the KINSOL memory block.infofp
– pointer to output file.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default value for
infofp
isstdout
.
Deprecated since version 6.2.0: Use
SUNLogger_SetInfoFilename()
instead.

int KINSetInfoHandlerFn(void *kin_mem, KINInfoHandlerFn ihfun, void *ih_data)
The function
KINSetInfoHandlerFn()
specifies the optional userdefined function to be used in handling informative (nonerror) messages. Arguments:
kin_mem
– pointer to the KINSOL memory block.ihfun
– is the user’s CC information handler function (see §8.4.6.3).ih_data
– pointer to user data passed toihfun
every time it is called.
 Return value:
KIN_SUCCESS
– The functionihfun
and data pointerih_data
have been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default internal information handler function directs informative (nonerror) messages to the file specified by the file pointer
infofp
(seeKINSetInfoFile()
above).

int KINSetPrintLevel(void *kin_mem, int printfl)
The function
KINSetPrintLevel()
specifies the level of verbosity of the output. Arguments:
kin_mem
– pointer to the KINSOL memory block.printfl
– flag indicating the level of verbosity. Must be one of:0 – no information is displayed.
1 – for each nonlinear iteration display the following information:
the scaled Euclidean \(\ell_2\) norm of the system function evaluated at the current iterate,
the scaled norm of the Newton step (only if using
KIN_NONE
), andthe number of function evaluations performed so far.
2 – display level 1 output and the following values for each iteration:
\(\F(u)\_{D_F}\) (only for
KIN_NONE
).\(\F(u)\_{D_F,\infty}\) (for
KIN_NONE
andKIN_LINESEARCH
).
3 – display level 2 output plus
additional values used by the global strategy (only if using
KIN_LINESEARCH
), andstatistical information for iterative linear solver modules.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentprintfl
had an illegal value.
 Notes:
The default value for
printfl
is \(0\).

int KINSetUserData(void *kin_mem, void *user_data)
The function
KINSetUserData()
specifies the pointer to userdefined memory that is to be passed to all usersupplied functions. Arguments:
kin_mem
– pointer to the KINSOL memory block.user_data
– pointer to the userdefined memory.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_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 toKINSetUserData()
must be made before the call to specify the linear solver module.

int KINSetNumMaxIters(void *kin_mem, long int mxiter)
The function
KINSetNumMaxIters()
specifies the maximum number of nonlinear iterations allowed. Arguments:
kin_mem
– pointer to the KINSOL memory block.mxiter
– maximum number of nonlinear iterations.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The maximum number of iterations was nonpositive.
 Notes:
The default value for
mxiter
isMXITER_DEFAULT
\(=200\).

int KINSetNoInitSetup(void *kin_mem, booleantype noInitSetup)
The function
KINSetNoInitSetup()
specifies whether an initial call to the preconditioner or Jacobian setup function should be made or not. Arguments:
kin_mem
– pointer to the KINSOL memory block.noInitSetup
– flag controlling whether an initial call to the preconditioner or Jacobian setup function is made (passSUNFALSE
) or not made (passSUNTRUE
).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default value for
noInitSetup
isSUNFALSE
, meaning that an initial call to the preconditioner or Jacobian setup function will be made. A call to this function is useful when solving a sequence of problems, in which the final preconditioner or Jacobian value from one problem is to be used initially for the next problem.

int KINSetNoResMon(void *kin_mem, booleantype noNNIResMon)
The function
KINSetNoResMon()
specifies whether or not the nonlinear residual monitoring scheme is used to control Jacobian updating Arguments:
kin_mem
– pointer to the KINSOL memory block.noNNIResMon
– flag controlling whether residual monitoring is used (passSUNFALSE
) or not used (passSUNTRUE
).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
When using a direct solver, the default value for
noNNIResMon
isSUNFALSE
, meaning that the nonlinear residual will be monitored.
Warning
Residual monitoring is only available for use with matrixbased linear solver modules.

int KINSetMaxSetupCalls(void *kin_mem, long int msbset)
The function
KINSetMaxSetupCalls()
specifies the maximum number of nonlinear iterations that can be performed between calls to the preconditioner or Jacobian setup function. Arguments:
kin_mem
– pointer to the KINSOL memory block.msbset
– maximum number of nonlinear iterations without a call to the preconditioner or Jacobian setup function. Pass 0 to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentmsbset
was negative.
 Notes:
The default value for
msbset
isMSBSET_DEFAULT=10
. The value ofmsbset
should be a multiple ofmsbsetsub
(seeKINSetMaxSubSetupCalls()
).

int KINSetMaxSubSetupCalls(void *kin_mem, long int msbsetsub)
The function
KINSetMaxSubSetupCalls()
specifies the maximum number of nonlinear iterations between checks by the residual monitoring algorithm. Arguments:
kin_mem
– pointer to the KINSOL memory block.msbsetsub
– maximum number of nonlinear iterations without checking the nonlinear residual. Pass 0 to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentmsbsetsub
was negative.
 Notes:
The default value for
msbsetsub
isMSBSET_SUB_DEFAULT
\(=5\). The value ofmsbset
(seeKINSetMaxSetupCalls()
) should be a multiple ofmsbsetsub
.
Warning
Residual monitoring is only available for use with matrixbased linear solver modules.

int KINSetEtaForm(void *kin_mem, int etachoice)
The function
KINSetEtaForm()
specifies the method for computing the value of the \(\eta\) coefficient used in the calculation of the linear solver convergence tolerance. Arguments:
kin_mem
– pointer to the KINSOL memory block.etachoice
– flag indicating the method for computing \(\eta\). The value must be one ofKIN_ETACHOICE1
,KIN_ETACHOICE2
, orKIN_ETACONSTANT
(see Chapter §8.2 for details).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentetachoice
had an illegal value.
 Notes:
The default value for
etachoice
isKIN_ETACHOICE1
. When using eitherKIN_ETACHOICE1
orKIN_ETACHOICE2
the safeguard\[\eta_n = \max(\eta_n, \eta_{\text{safe}})\]is applied when \(\eta_{\text{safe}} > 0.1\). For
KIN_ETACHOICE1
\[\eta_{\text{safe}} = \eta_{n1}^{\frac{1+\sqrt{5}}{2}}\]and for
KIN_ETACHOICE2
\[\eta_{\text{safe}} = \gamma \eta_{n1}^\alpha\]where \(\gamma\) and \(\alpha\) can be set with
KINSetEtaParams()
.The following safeguards are always applied when using either
KIN_ETACHOICE1
orKIN_ETACHOICE2
so that \(\eta_{\text{min}} \leq \eta_n \leq\eta_{\text{max}}\):\[\begin{split}\begin{aligned} \eta_n &= \max(\eta_n, \eta_{\text{min}}) \\ \eta_n &= \min(\eta_n, \eta_{\text{max}}) \end{aligned}\end{split}\]where \(\eta_{\text{min}} = 10^{4}\) and \(\eta_{\text{max}} = 0.9\).

int KINSetEtaConstValue(void *kin_mem, realtype eta)
The function
KINSetEtaConstValue()
specifies the constant value for \(\eta\) in the caseetachoice = KIN_ETACONSTANT
. Arguments:
kin_mem
– pointer to the KINSOL memory block.eta
– constant value for \(\eta\). Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumenteta
had an illegal value
 Notes:
The default value for
eta
is \(0.1\). The legal values are \(0.0 <\)eta
\(\le 1.0\).

int KINSetEtaParams(void *kin_mem, realtype egamma, realtype ealpha)
The function
KINSetEtaParams()
specifies the parameters \(\gamma\) and \(\alpha\) in the formula for \(\eta\), in the caseetachoice = KIN_ETACHOICE2
. Arguments:
kin_mem
– pointer to the KINSOL memory block.egamma
– value of the \(\gamma\) parameter. Pass \(0.0\) to indicate the default.ealpha
– value of the \(\alpha\) parameter. Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional values have been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– One of the argumentsegamma
orealpha
had an illegal value.
 Notes:
The default values for
egamma
andealpha
are \(0.9\) and \(2.0\), respectively. The legal values are \(0.0 <\)egamma
\(\le 1.0\) and \(1.0<\)ealpha
\(\le 2.0\).

int KINSetResMonConstValue(void *kin_mem, realtype omegaconst)
The function
KINSetResMonConstValue()
specifies the constant value for \(\omega\) when using residual monitoring. Arguments:
kin_mem
– pointer to the KINSOL memory block.omegaconst
– constant value for \(\omega\). Passing \(0.0\) results in using Eqn. (8.4).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentomegaconst
had an illegal value
 Notes:
The default value for
omegaconst
is \(0.9\). The legal values are \(0.0 <\)omegaconst
\(< 1.0\).

int KINSetResMonParams(void *kin_mem, realtype omegamin, realtype omegamax)
The function
KINSetResMonParams()
specifies the parameters \(\omega_{min}\) and \(\omega_{max}\) in the formula (8.4) for \(\omega\). Arguments:
kin_mem
– pointer to the KINSOL memory block.omegamin
– value of the \(\omega_{min}\) parameter. Pass \(0.0\) to indicate the default.omegamax
– value of the \(\omega_{max}\) parameter. Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional values have been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– One of the argumentsomegamin
oromegamax
had an illegal value.
 Notes:
The default values for
omegamin
andomegamax
are \(0.00001\) and \(0.9\), respectively. The legal values are \(0.0 <\)omegamin
\(<\)omegamax
\(< 1.0\).
Warning
Residual monitoring is only available for use with matrixbased linear solver modules.

int KINSetNoMinEps(void *kin_mem, booleantype noMinEps)
The function
KINSetNoMinEps()
specifies a flag that controls whether or not the value of \(\epsilon\), the scaled linear residual tolerance, is bounded from below. Arguments:
kin_mem
– pointer to the KINSOL memory block.noMinEps
– flag controlling the bound on \(\epsilon\). IfSUNFALSE
is passed the value of \(\epsilon\) is constrained and ifSUNTRUE
is passed then \(\epsilon\) is not constrained.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default value for
noMinEps
isSUNFALSE
, meaning that a positive minimum value, equal to \(0.01\), is applied to \(\epsilon\) (seeKINSetFuncNormTol()
below).

int KINSetMaxNewtonStep(void *kin_mem, realtype mxnewtstep)
The function
KINSetMaxNewtonStep()
specifies the maximum allowable scaled length of the Newton step. Arguments:
kin_mem
– pointer to the KINSOL memory block.mxnewtstep
– maximum scaled step length \((\geq 0.0)\). Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The input value was negative.
 Notes:
The default value of
mxnewtstep
is \(1000\, \ u_0 \_{D_u}\), where \(u_0\) is the initial guess.

int KINSetMaxBetaFails(void *kin_mem, realtype mxnbcf)
The function
KINSetMaxBetaFails()
specifies the maximum number of \(\beta\)condition failures in the linesearch algorithm. Arguments:
kin_mem
– pointer to the KINSOL memory block.mxnbcf
– maximum number of \(\beta\) condition failures. Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
–mxnbcf
was negative.
 Notes:
The default value of
mxnbcf
isMXNBCF_DEFAULT
\(=10\).

int KINSetRelErrFunc(void *kin_mem, realtype relfunc)
The function
KINSetRelErrFunc()
specifies the relative error in computing \(F(u)\), which is used in the difference quotient approximation to the Jacobian matrix [see Eq. (8.6) ] or the Jacobianvector product [see Eq. (8.8) ]. The value stored is \(\sqrt{\texttt{relfunc}}\). Arguments:
kin_mem
– pointer to the KINSOL memory block.relfunc
– relative error in \(F(u)\) (\(\texttt{relfunc} \geq 0.0\)). Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The relative error was negative.
 Notes:
The default value for
relfunc
is \(U\) = unit roundoff.

int KINSetFuncNormTol(void *kin_mem, realtype fnormtol)
The function
KINSetFuncNormTol()
specifies the scalar used as a stopping tolerance on the scaled maximum norm of the system function \(F(u)\). Arguments:
kin_mem
– pointer to the KINSOL memory block.fnormtol
– tolerance for stopping based on scaled function norm \((\geq 0.0)\). Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The tolerance was negative.
 Notes:
The default value for
fnormtol
is (unit roundoff) \(^{1/3}\).

int KINSetScaledStepTol(void *kin_mem, realtype scsteptol)
The function
KINSetScaledStepTol()
specifies the scalar used as a stopping tolerance on the minimum scaled step length. Arguments:
kin_mem
– pointer to the KINSOL memory block.scsteptol
– tolerance for stopping based on scaled step length \((\geq 0.0)\). Pass \(0.0\) to indicate the default.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The tolerance was nonpositive.
 Notes:
The default value for
scsteptol
is (unit roundoff) \(^{2/3}\).

int KINSetConstraints(void *kin_mem, N_Vector constraints)
The function
KINSetConstraints()
specifies a vector that defines inequality constraints for each component of the solution vector \(u\). Arguments:
kin_mem
– pointer to the KINSOL memory block.constraints
– vector of constraint flags. Ifconstraints[i]
is\(0.0\) then no constraint is imposed on \(u_i\).
\(1.0\) then \(u_i\) will be constrained to be \(u_i \ge 0.0\).
\(1.0\) then \(u_i\) will be constrained to be \(u_i \le 0.0\).
\(2.0\) then \(u_i\) will be constrained to be \(u_i > 0.0\).
\(2.0\) then \(u_i\) will be constrained to be \(u_i < 0.0\).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The constraint 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. If aNULL
vector is supplied, constraint checking will be disabled. The function creates a private copy of the constraints vector. Consequently, the usersupplied vector can be freed after the function call, and the constraints can only be changed by calling this function.

int KINSetSysFunc(void *kin_mem, KINSysFn func)
The function
KINSetSysFunc()
specifies the userprovided function that evaluates the nonlinear system function \(F(u)\) or \(G(u)\). Arguments:
kin_mem
– pointer to the KINSOL memory block.func
– usersupplied function that evaluates \(F(u)\) (or \(G(u)\) for fixedpoint iteration).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentfunc
wasNULL
.
 Notes:
The nonlinear system function is initially specified through
KINInit()
. The option of changing the system function is provided for a user who wishes to solve several problems of the same size but with different functions.

int KINSetReturnNewest(void *kin_mem, booleantype ret_newest)
The function
KINSetReturnNewest()
specifies if the fixed point iteration should return the newest iteration or the iteration consistent with the last function evaluation. Arguments:
kin_mem
– pointer to the KINSOL memory block.ret_newest
–SUNTRUE
– return the newest iteration.SUNFALSE
– return the iteration consistent with the last function evaluation.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
The default value of
ret_newest
isSUNFALSE
.

int KINSetDamping(void *kin_mem, realtype beta)
The function
KINSetDamping()
specifies the value of the damping parameter in the fixed point or Picard iteration. Arguments:
kin_mem
– pointer to the KINSOL memory block.beta
– the damping parameter value \(0 < beta \leq 1.0\).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentbeta
was zero or negative.
 Notes:
This function sets the damping parameter value, which needs to be greater than zero and less than one if damping is to be used. A value \(\geq 1\) disables damping. The default value of
beta
is 1.0, indicating no damping. To set the damping parameter used in Anderson acceleration seeKINSetDampingAA()
. With the fixed point iteration the difference between successive iterations is used to determine convergence. As such, when damping is enabled, the tolerance used to stop the fixed point iteration is scaled bybeta
to account for the effects of damping. Ifbeta
is extremely small (close to zero), this can lead to an excessively tight tolerance.

int KINSetMAA(void *kin_mem, long int maa)
The function
KINSetMAA()
specifies the size of the subspace used with Anderson acceleration in conjunction with Picard or fixedpoint iteration. Arguments:
kin_mem
– pointer to the KINSOL memory block.maa
– subspace size for various methods. A value of 0 means no acceleration, while a positive value means acceleration will be done.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentmaa
was negative.
 Notes:
This function sets the subspace size, which needs to be \(> 0\) if Anderson Acceleration is to be used. It also allocates additional memory necessary for Anderson Acceleration. The default value of
maa
is 0, indicating no acceleration. The value ofmaa
should always be less thanmxiter
. This function MUST be called before callingKINInit()
. If the user calls the function KINSetNumMaxIters, that call should be made before the call to KINSetMAA, as the latter uses the value ofmxiter
.

int KINSetDampingAA(void *kin_mem, realtype beta)
The function
KINSetDampingAA()
specifies the value of the Anderson acceleration damping paramter. Arguments:
kin_mem
– pointer to the KINSOL memory block.beta
– the damping parameter value \(0 < beta \leq 1.0\).
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentbeta
was zero or negative.
 Notes:
This function sets the damping parameter value, which needs to be greater than zero and less than one if damping is to be used. A value \(\geq 1\) disables damping. The default value of
beta
is 1.0, indicating no damping. When delaying the start of Anderson acceleration withKINSetDelayAA()
, useKINSetDamping()
to set the damping parameter in the fixed point or Picard iterations before Anderson acceleration begins. When using Anderson acceleration without delay, the value provided toKINSetDampingAA()
is applied to all iterations and any value provided toKINSetDamping()
is ignored.

int KINSetDelayAA(void *kin_mem, long int delay)
The function
KINSetDelayAA()
specifies the number of iterations to delay the start of Anderson acceleration. Arguments:
kin_mem
– pointer to the KINSOL memory block.delay
– the number of iterations to delay Anderson acceleration.
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentdelay
was less than zero.
 Notes:
The default value of
delay
is 0, indicating no delay.

int KINSetOrthAA(void *kin_mem, int orthaa)
The function
KINSetOrthAA()
specifies the orthogonalization routine to be used in the QR factorization portion of Anderson acceleration. Arguments:
kin_mem
– pointer to the KINSOL memory block.orthaa
– the orthogonalization routine parameter. Can be set to any ofthe following
KIN_ORTH_MGS
– Modified Gram Schmidt (default)KIN_ORTH_ICWY
– Inverse Compact WY Modified Gram SchmidtKIN_ORTH_CGS2
– Classical Gram Schmidt with Reorthogonalization (CGS2)KIN_ORTH_DCGS2
– Classical Gram Schmidt with Delayed Reorthogonlization
 Return value:
KIN_SUCCESS
– The optional value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_ILL_INPUT
– The argumentorthaa
was not one of the predefined orthogonalization routines defined in KINSOL.
Note
This function must be called before calling
KINInit()
.An example of how to use this function can be found in
examples/kinsol/serial/kinAnalytic_fp.c
8.4.5.4.1. Linear solver interface optional input functions
For matrixbased linear solver modules, the KINLS solver interface needs a
function to compute an approximation to the Jacobian matrix \(J(u)\). This
function must be of type KINLsJacFn
. The user can supply a Jacobian
function, or if using the SUNMATRIX_DENSE or
SUNMATRIX_BAND modules for \(J\) can use the default
internal difference quotient approximation that comes with the KINLS solver. To
specify a usersupplied Jacobian function jac
, KINLS provides the function
KINSetJacFn()
. The KINLS 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 KINSetUserData()
.

int KINSetJacFn(void *kin_mem, KINLsJacFn jac)
The function
KINSetJacFn()
specifies the Jacobian approximation function to be used for a matrixbased solver within the KINLS interface. Arguments:
kin_mem
– pointer to the KINSOL solver object.jac
– userdefined Jacobian approximation function. SeeKINLsJacFn
for more details.
 Return value:
KINLS_SUCCESS
– The optional value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver interface has not been initialized.
 Notes:
This function must be called after the KINLS linear solver interface has been initialized through a call to
KINSetLinearSolver()
. By default, KINLS 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
KINDlsSetJacFn()
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 matrixfree linear solver modules, the KINLS linear solver interface requires a function to compute an approximation to the product between the Jacobian matrix \(J(u)\) and a vector \(v\). The user can supply his/her own Jacobiantimesvector approximation function, or use the internal difference quotient approximation that comes with the KINLS solver interface.
A userdefined Jacobianvector function must be of type KINLsJacTimesVecFn
and can be specified through a call to KINLsSetJacTimesVecFn()
(see
§8.4.6.5 for specification details). The pointer
user_data
received through KINSetUserData()
(or a pointer to NULL
if
user_data
was not specified) is passed to the Jacobiantimesvector function
jtimes
each time it is 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 KINSetJacTimesVecFn(void *kin_mem, KINLsJacTimesVecFn jtimes)
The function
KINSetJacTimesVecFn()
specifies the Jacobianvector product function. Arguments:
kin_mem
– pointer to the KINSOL memory block.jtimes
– userdefined Jacobianvector product function.
 Return value:
KINLS_SUCCESS
– The optional value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.KINLS_SUNLS_FAIL
– An error occurred when setting up the system matrixtimesvector routines in the SUNLINSOL object used by the KINLS interface.
 Notes:
The default is to use an internal difference quotient for
jtimes
. IfNULL
is passed asjtimes
, this default is used. This function must be called after the KINLS linear solver interface has been initialized through a call toKINSetLinearSolver()
. The function typeKINLsJacTimesVecFn
is described in §8.4.6.5. The previous routineKINSpilsSetJacTimesVecFn()
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 internal difference quotient the user may optionally supply an
alternative system function for use in the Jacobianvector product approximation
by calling KINSetJacTimesVecSysFn()
. The alternative system function
should compute a suitable (and differentiable) approximation of the system
function provided to KINInit()
. For example, as done in
[44] when solving the nonlinear systems that arise in the
implicit integration of ordinary differential equations, the alternative
function may use lagged values when evaluating a nonlinearity to avoid
differencing a potentially nondifferentiable factor.

int KINSetJacTimesVecSysFn(void *kin_mem, KINSysFn jtimesSysFn)
The function
KINSetJacTimesVecSysFn()
specifies an alternative system function for use in the internal Jacobianvector product difference quotient approximation. Arguments:
kin_mem
– pointer to the KINSOL memory block.jtimesSysFn
– is the CC function which computes the alternative system function to use in Jacobianvector product difference quotient approximations. This function has the formfunc(u, fval, user_data)
. (For full details see §8.4.6.1.)
 Return value:
KINLS_SUCCESS
– The optional value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.KINLS_ILL_INPUT
– The internal difference quotient approximation is disabled.
 Notes:
The default is to use the system function provided to
KINInit()
in the internal difference quotient. If the input system function isNULL
, the default is used. This function must be called after the KINLS linear solver interface has been initialized through a call toKINSetLinearSolver()
.
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 KINLS
using the function KINSetPreconditioner()
. 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 §8.4.6.
The user data pointer received through
KINSetUserData()
(or a pointer to NULL
if user data 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 KINSetPreconditioner(void *kin_mem, KINLsPrecSetupFn psetup, KINLsPrecSolveFn psolve)
The function
KINSetPreconditioner()
specifies the preconditioner setup and solve functions. Arguments:
kin_mem
– pointer to the KINSOL solver object.psetup
– userdefined function to set up the preconditioner. SeeKINLsPrecSetupFn
for more details. PassNULL
if no setup is necessary.psolve
– userdefined preconditioner solve function. SeeKINLsPrecSolveFn
for more details.
 Return value:
KINLS_SUCCESS
– The optional values have been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.KINLS_SUNLS_FAIL
– An error occurred when setting up preconditioning in theSUNLinearSolver
object used by the KINLS interface.
 Notes:
The default is
NULL
for both arguments (i.e., no preconditioning). This function must be called after the KINLS linear solver interface has been initialized through a call toKINSetLinearSolver()
.
Warning
The previous routine
KINSpilsSetPreconditioner()
is now a wrapper for this routine, and may still be used for backwardcompatibility. However, this will be removed in future releases, so we recommend that users transition to the new routine name soon.
8.4.5.5. Optional output functions
KINSOL provides an extensive list of functions that can be used to obtain solver
performance information. Table 8.2
lists all optional output functions in KINSOL, which are then described in
detail in the remainder of this section, beginning with those for the main
KINSOL solver and continuing with those for the KINLS linear solver interface.
Where the name of an output from a linear solver module would otherwise conflict
with the name of an optional output from the main solver, a suffix LS
(for
Linear Solver) has been added here (e.g., lenrwLS
).
Optional output 
Function name 

KINSOL main solver 

Size of KINSOL real and integer workspaces 

Number of function evaluations 

Number of nonlinear iterations 

Number of \(\beta\)condition failures 

Number of backtrack operations 

Scaled norm of \(F\) 

Scaled norm of the step 

User data pointer 

Print all statistics 

Name of constant associated with a return flag 

KINLS linear solver interface 

Stored Jacobian of the nonlinear system 

Nonlinear iteration number at which the Jacobian was evaluated 

Size of real and integer workspaces 

No. of Jacobian evaluations 

No. of \(F\) calls for D.Q. Jacobian[vector] evals. 

No. of linear iterations 

No. of linear convergence failures 

No. of preconditioner evaluations 

No. of preconditioner solves 

No. of Jacobianvector product evaluations 

Last return from a KINLS function 

Name of constant associated with a return flag 
8.4.5.5.1. Main solver optional output functions
KINSOL 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 and solver performance statistics. These optional output functions are described next.

int KINGetWorkSpace(void *kin_mem, long int lenrw, long int leniw)
The function
KINGetWorkSpace()
returns the KINSOL integer and real workspace sizes. Arguments:
kin_mem
– pointer to the KINSOL memory block.lenrw
– the number ofrealtype
values in the KINSOL workspace.leniw
– the number of integer values in the KINSOL workspace.
 Return value:
KIN_SUCCESS
– The optional output values have been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
 Notes:
KINSOL solver In terms of the problem size \(N\), the actual size of the real workspace is \(17 + 5 N\)
realtype
words. The real workspace is increased by an additional \(N\) words if constraint checking is enabled (seeKINSetConstraints()
).The actual size of the integer workspace (without distinction between
int
andlong int
) is \(22 + 5 N\) (increased by \(N\) if constraint checking is enabled).

int KINGetNumFuncEvals(void *kin_mem, long int nfevals)
The function
KINGetNumFuncEvals()
returns the number of evaluations of the system function. Arguments:
kin_mem
– pointer to the KINSOL memory block.nfevals
– number of calls to the usersupplied function that evaluates \(F(u)\).
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetNumNonlinSolvIters(void *kin_mem, long int nniters)
The function
KINGetNumNonlinSolvIters()
returns the number of nonlinear iterations. Arguments:
kin_mem
– pointer to the KINSOL memory block.nniters
– number of nonlinear iterations.
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetNumBetaCondFails(void *kin_mem, long int nbcfails)
The function
KINGetNumBetaCondFails()
returns the number of \(\beta\)condition failures. Arguments:
kin_mem
– pointer to the KINSOL memory block.nbcfails
– number of \(\beta\) condition failures.
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetNumBacktrackOps(void *kin_mem, long int nbacktr)
The function
KINGetNumBacktrackOps()
returns the number of backtrack operations (step length adjustments) performed by the line search algorithm. Arguments:
kin_mem
– pointer to the KINSOL memory block.nbacktr
– number of backtrack operations.
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetFuncNorm(void *kin_mem, realtype fnorm)
The function
KINGetFuncNorm()
returns the scaled Euclidean \(\ell_2\) norm of the nonlinear system function \(F(u)\) evaluated at the current iterate. Arguments:
kin_mem
– pointer to the KINSOL memory block.fnorm
– current scaled norm of \(F(u)\).
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetStepLength(void *kin_mem, realtype steplength)
The function
KINGetStepLength()
returns the scaled Euclidean \(\ell_2\) norm of the step used during the previous iteration. Arguments:
kin_mem
– pointer to the KINSOL memory block.steplength
– scaled norm of the Newton step.
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.

int KINGetUserData(void *kin_mem, void **user_data)
The function
KINGetUserData()
returns the user data pointer provided toKINSetUserData()
. Arguments:
kin_mem
– pointer to the KINSOL memory block.user_data
– memory reference to a user data pointer.
 Return value:
KIN_SUCCESS
– The optional output value has been successfully set.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.
New in version 6.3.0.

int KINPrintAllStats(void *cvode_mem, FILE *outfile, SUNOutputFormat fmt)
The function
KINPrintAllStats()
outputs all of the nonlinear solver, linear solver, and other statistics. Arguments:
kin_mem
– pointer to the KINSOL 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:
KIN_SUCCESS
– The output was successfully.KIN_MEM_NULL
– Thekin_mem
pointer isNULL
.KIN_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 *KINGetReturnFlagName(int flag)
The function
KINGetReturnFlagName()
returns the name of the KINSOL constant corresponding toflag
. Arguments:
flag
– return flag from a KINSOL function.
 Return value:
A string containing the name of the corresponding constant
8.4.5.5.2. KINLS linear solver interface optional output functions
The following optional outputs are available from the KINLS modules:

int KINGetJac(void *kin_mem, SUNMatrix *J)
Returns the internally stored copy of the Jacobian matrix of the nonlinear system function.
 Parameters
kin_mem – the KINSOL solver object
J – the Jacobian matrix
 Return values
KINLS_SUCCESS – the output value has been successfully set
KINLS_MEM_NULL –
kin_mem
wasNULL
KINLS_LMEM_NULL – the linear solver interface has not been initialized
Warning
With linear solvers that overwrite the input Jacobian matrix as part of the linear solver setup (e.g., performing an inplace LU factorization) the matrix returned by
KINGetJac()
may differ from the matrix returned by the last Jacobian evaluation.Warning
This function is provided for debugging purposes and the values in the returned matrix should not be altered.

int KINGetJacNumIters(void *kin_mem, sunrealtype *nni_J)
Returns the nonlinear iteration number at which the Jacobian was evaluated.
 Parameters
kin_mem – the KINSOL memory structure
nni_J – the nonlinear iteration number
 Return values
KINLS_SUCCESS – the output value has been successfully set
KINLS_MEM_NULL –
kin_mem
wasNULL
KINLS_LMEM_NULL – the linear solver interface has not been initialized

int KINGetLinWorkSpace(void *kin_mem, long int *lenrwLS, long int *leniwLS)
The function
KINGetLinWorkSpace()
returns the sizes of the real and integer workspaces used by the KINLS linear solver interface. Arguments:
kin_mem
– pointer to the KINSOL solver object.lenrwLS
– the number of real values in the KINLS workspace.leniwLS
– the number of integer values in the KINLS workspace.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS 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 KINLS is not included in this report.
Warning
The previous routines
KINDlsGetWorkspace()
andKINSpilsGetWorkspace()
are now deprecated.

int KINGetNumJacEvals(void *kin_mem, long int *njevals)
The function
KINGetNumJacEvals()
returns the cumulative number of calls to the KINLS Jacobian approximation function. Arguments:
kin_mem
– pointer to the KINSOL solver object.njevals
– the cumulative number of calls to the Jacobian function total so far.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINDlsGetNumJacEvals()
is now deprecated,

int KINGetNumLinFuncEvals(void *kin_mem, long int *nrevalsLS)
The function
KINGetNumLinResEvals()
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:
kin_mem
– pointer to the KINSOL solver object.nrevalsLS
– the cumulative number of calls to the user residual function.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS 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
KINDlsGetNumRhsEvals()
andKINSpilsGetNumRhsEvals()
are now deprecated.

int KINGetNumLinIters(void *kin_mem, long int *nliters)
The function
KINGetNumLinIters()
returns the cumulative number of linear iterations. Arguments:
kin_mem
– pointer to the KINSOL solver object.nliters
– the current number of linear iterations.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINSpilsGetNumLinIters()
is now deprecated.

int KINGetNumLinConvFails(void *kin_mem, long int *nlcfails)
The function
KINGetNumLinConvFails()
returns the cumulative number of linear convergence failures. Arguments:
kin_mem
– pointer to the KINSOL solver object.nlcfails
– the current number of linear convergence failures.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINSpilsGetNumConvFails()
is now deprecated.

int KINGetNumPrecEvals(void *kin_mem, long int *npevals)
The function
KINGetNumPrecEvals()
returns the cumulative number of preconditioner evaluations, i.e., the number of calls made topsetup
. Arguments:
kin_mem
– pointer to the KINSOL solver object.npevals
– the cumulative number of calls topsetup
.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINSpilsGetNumPrecEvals()
is now deprecated.

int KINGetNumPrecSolves(void *kin_mem, long int *npsolves)
The function
KINGetNumPrecSolves()
returns the cumulative number of calls made to the preconditioner solve function,psolve
. Arguments:
kin_mem
– pointer to the KINSOL solver object.npsolves
– the cumulative number of calls topsolve
.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINSpilsGetNumPrecSolves()
is now deprecated.

int KINGetNumJtimesEvals(void *kin_mem, long int *njvevals)
The function
KINGetNumJtimesEvals()
returns the cumulative number of calls made to the Jacobianvector product function,jtimes
. Arguments:
kin_mem
– pointer to the KINSOL solver object.njvevals
– the cumulative number of calls tojtimes
.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
Warning
The previous routine
KINSpilsGetNumJtimesEvals()
is now deprecated.

int KINGetLastLinFlag(void *kin_mem, long int *lsflag)
The function
KINGetLastLinFlag()
returns the last return value from an KINLS routine. Arguments:
kin_mem
– pointer to the KINSOL solver object.lsflag
– the value of the last return flag from an KINLS function.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer isNULL
.KINLS_LMEM_NULL
– The KINLS linear solver has not been initialized.
 Notes:
If the KINLS setup function failed (i.e.,
KINSolve()
returnedKIN_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 KINLS setup function failed when using another
SUNLinearSolver
object, thenlsflag
will beSUNLS_PSET_FAIL_UNREC
,SUNLS_ASET_FAIL_UNREC
, orSUNLS_PACKAGE_FAIL_UNREC
.If the KINLS solve function failed (
KINSolve()
returnedKIN_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
KINDlsGetLastFlag()
andKINSpilsGetLastFlag()
are now deprecated.

char *KINGetLinReturnFlagName(long int lsflag)
The function
KINGetLinReturnFlagName()
returns the name of the KINLS constant corresponding tolsflag
. Arguments:
flag
– the flag returned by a call to an KINSOL 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
KINDlsGetReturnFlagName()
andKINSpilsGetReturnFlagName()
are now deprecated.
8.4.6. Usersupplied functions
The usersupplied functions consist of one function defining the nonlinear system, (optionally) a function that handles error and warning messages, (optionally) a function that handles informational messages, (optionally) one or two functions that provides Jacobianrelated information for the linear solver, and (optionally) one or two functions that define the preconditioner for use in any of the Krylov iterative algorithms.
8.4.6.1. Problem defining function
The user must provide a function of type KINSysFn
defined as follows:

typedef int (*KINSysFn)(N_Vector u, N_Vector fval, void *user_data)
This function computes the \(F(u)\) (or \(G(u)\) for fixedpoint iteration and Anderson acceleration) for a given value of the vector \(u\).
 Arguments:
u
– is the current value of the dependent variable vector, \(u\)fval
– is the output vector \(F(u)\)user_data
– is a pointer to user data, the same as theuser_data
pointer parameter passed toKINSetUserData()
 Return value:
An
KINSysFn
function type should return a value of \(0\) if successful, a positive value if a recoverable error occurred (in which case KINSOL will attempt to correct), 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:
Allocation of memory for
fval
is handled within KINSOL.
8.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 KINSetErrFile()
), the
user may provide a function of type KINErrHandlerFn
to process any
such messages. The function type KINErrHandlerFn
is defined as
follows:

typedef void (*KINErrHandlerFn)(int error_code, const char *module, const char *function, char *msg, void *user_data)
This function processes error and warning messages from KINSOL and its submodules.
 Arguments:
error_code
– is the error codemodule
– is the name of the KINSOL module reporting the errorfunction
– is the name of the function in which the error occurredeH_data
– is a pointer to user data, the same as theeh_data
parameter passed toKINSetErrHandlerFn()
 Return value:
This function has no return value.
 Notes:
error_code
is negative for errors and positive (KIN_WARNING
) for warnings. If a function that returns a pointer to memory encounters an error, it setserror_code
to 0.
8.4.6.3. Informational message handler function
As an alternative to the default behavior of directing informational (meaning nonerror) messages
to the file pointed to by infofp
(see KINSetInfoFile()
), the user may
provide a function of type KINInfoHandlerFn
to process any such messages.
The function type KINInfoHandlerFn
is defined as follows:

typedef void (*KINInfoHandlerFn)(const char *module, const char *function, char *msg, void *ih_data)
This function processes error and warning messages from KINSOL and its submodules.
 Arguments:
error_code
– is the error codemodule
– is the name of the KINSOL module reporting the errorfunction
– is the name of the function in which the error occurredih_data
– is a pointer to user data, the same as theih_data
parameter passed toKINSetInfoHandlerFn()
 Return value:
This function has no return value.
8.4.6.4. Jacobian construction (matrixbased linear solvers)
If a matrixbased linear solver module is used (i.e. a nonNULL
SUNMatrix
object was supplied to KINSetLinearSolver()
), the user may
provide a function of type KINLsJacFn
defined as follows:

typedef int (*KINLsJacFn)(N_Vector u, N_Vector fu, SUNMatrix J, void *user_data, N_Vector tmp1, N_Vector tmp2)
This function computes the Jacobian matrix \(J(u)\) (or an approximation to it).
 Arguments:
u
– is the current (unscaled) iterate.fu
– is the current value of the vector, \(F(u)\).J
– is the output (approximate) Jacobian matrix (of typeSUNMatrix
), \(F'(u)\).user_data
 is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.tmp1
,tmp2
, – are pointers to memory allocated for variables of typeN_Vector
which can be used byKINLsJacFn
function as temporary storage or work space.
 Return value:
An
KINLsJacFn
should return \(0\) if successful, or a nonzero value otherwise. 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(u)\) is zeroed out prior to calling the usersupplied Jacobian function so only nonzero elements need to be loaded intoJ
.If the user’s
KINLsJacFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These quantities may include the scale vectors and the unit roundoff. To obtain the scale vectors, the user will need to add touser_data
pointers tou_scale
and/orf_scale
as needed. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.dense:
A usersupplied dense Jacobian function must load the
N
\(\times\)N
dense matrixJ
with an approximation to the Jacobian matrix \(J(u)\) at the point (u
). 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 matrixJ
(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 ofJ
(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
N
\(\times\)N
banded matrixJ
with an approximation to the Jacobian matrix \(J(u)\) at the point (u
). 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 matrixJ
, 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 ofJ
, 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
N
\(\times\)N
compressedsparsecolumn or compressedsparserow matrixJ
with an approximation to the Jacobian matrix \(J(u)\) at the point (u
). Storage forJ
already exists on entry to this function, although the user should ensure that sufficient space is allocated inJ
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
KINDlsJacFn()
is identical toKINLsJacFn
, 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.
8.4.6.5. Jacobianvector product (matrixfree linear solvers)
If a matrixfree linear solver is to be used (i.e., a NULL
valued
SUNMatrix
was supplied to KINSetLinearSolver()
), the user may
provide a function of type KINLsJacTimesVecFn
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 (*KINLsJacTimesVecFn)(N_Vector v, N_Vector Jv, N_Vector u, booleantype *new_u, void *user_data)
This function computes the product \(J v\) (or an approximation to it).
 Arguments:
v
– is the vector by which the Jacobian must be multplied to the right.Jv
– is the computed output vector.u
– is the current value of the dependent variable vector.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.
 Return value:
The value returned by the Jacobiantimesvector function should be 0 if successful. If a recoverable failure occurred, the return value should be positive. In this case, KINSOL will attempt to correct by calling the preconditioner setup function. If this information is current, KINSOL halts. If the Jacobiantimesvector function encounters an unrecoverable error, it should return a negative value, prompting KINSOL to halt.
 Notes:
If a userdefined routine is not given, then an internal
jtimes
function, using a difference quotient approximation, is used.This function must return a value of \(J*v\) that uses the current value of \(J\), i.e. as evaluated at the current \(u\).
If the user’s
KINLsJacTimesVecFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These might include the scale vectors and the unit roundoff. To obtain the scale vectors, the user will need to add touser_data
pointers tou_scale
and/orf_scale
as needed. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.
Warning
The previous function type
KINSpilsJacTimesVecFn
is identical toKINLsJacTimesVecFn
, and may still be used for backwardcompatibility. However, this will be removed in future releases, so we recommend that users transition to the new function type name soon.
8.4.6.6. 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 the preconditioner matrix which approximates
(at least crudely) the Jacobian matrix \(J = F'(u)\). This function must be
of type KINLsPrecSolveFn
, defined as follows:

typedef int (*KINLsPrecSolveFn)(N_Vector u, N_Vector uscale, N_Vector fval, N_Vector fscale, N_Vector v, void *user_data)
This function solves the preconditioning system \(Pz = r\).
 Arguments:
u
– is the current (unscaled) value of the iterate.uscale
– is a vector containing diagonal elements of the scaling matrixu
fval
– is the vector \(F(u)\) evaluated atu
fscale
– is a vector containing diagonal elements of the scaling matrix forfval
v
– on inpuut,v
is set to the righthand side vector of the linear system,r
. On output,v
must contain the solutionz
of the linear system \(Pz=r\)user_data
– is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.
 Return value:
The value returned by the preconditioner solve function should be 0 if successful, positive for a recoverable error, or negative for an unrecoverable error.
 Notes:
If the preconditioner solve function fails recoverably and if the preconditioner information (set by the preconditioner setup function) is out of date, KINSOL attempts to correct by calling the setup function. If the preconditioner data is current, KINSOL halts.
8.4.6.7. 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 KINLsPrecSetupFn
, defined as follows:

typedef int (*KINLsPrecSetupFn)(N_Vector u, N_Vector uscale, N_Vector fval, N_Vector fscale, void *user_data)
This function evaluates and/or preprocesses Jacobianrelated data needed by the preconditioner solve function.
 Arguments:
u
– is the current (unscaled) value of the iterate.uscale
– is a vector containing diagonal elements of the scaling matrixu
fval
– is the vector \(F(u)\) evaluated atu
fscale
– is a vector containing diagonal elements of the scaling matrix forfval
user_data
– is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.
 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:
The usersupplied preconditioner setup subroutine should compute the right preconditioner matrix \(P\) (stored in the memory block referenced by the
user_data
pointer) used to form the scaled preconditioned linear system\[(D_F J(u) P^{1} D_u^{1}) (D_u P x) =  D_F F(u) \, ,\]where \(D_u\) and \(D_F\) denote the diagonal scaling matrices whose diagonal elements are stored in the vectors
uscale
andfscale
, respectively.The preconditioner setup routine will not be called prior to every call made to the preconditioner solve function, but will instead be called only as often as necessary to achieve convergence of the Newton iteration.
If the user’s
KINLsPrecSetupFn
function uses difference quotient approximations, it may need to access quantities not in the call list. These might include the scale vectors and the unit roundoff. To obtain the scale vectors, the user will need to add touser_data
pointers tou_scale
and/orf_scale
as needed. The unit roundoff can be accessed asUNIT_ROUNDOFF
defined insundials_types.h
.If the preconditioner solve routine requires no preparation, then a preconditioner setup function need not be given.
8.4.7. A parallel bandblockdiagonal preconditioner module
The efficiency of Krylov iterative methods for the solution of linear systems
can be greatly enhanced through preconditioning. For problems in which the user
cannot define a more effective, problemspecific preconditioner, KINSOL provides
a bandblockdiagonal preconditioner module KINBBDPRE, to be used with the
parallel N_Vector
module described in §9.10.
This module provides a preconditioner matrix for KINSOL that is blockdiagonal
with banded blocks. The blocking corresponds to the distribution of the
dependent variable vector \(u\) amongst the processes. Each preconditioner
block is generated from the Jacobian of the local part (associated with the
current process) of a given function \(G(u)\) approximating \(F(u)\)
(\(G = F\) is allowed). The blocks are generated by each process via a
difference quotient scheme, utilizing a specified band structure. This structure
is given by upper and lower halfbandwidths, mudq
and mldq
, defined as
the number of nonzero diagonals above and below the main diagonal,
respectively. However, from the resulting approximate Jacobain blocks, only a
matrix of bandwidth mukeep
\(+\) mlkeep
\(+ 1\) is retained.
Neither pair of parameters need be the true halfbandwidths of the Jacobian 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 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
to see whether the lower cost of narrower band matrices offsets the loss of
accuracy in the blocks.
The KINBBDPRE module calls two userprovided functions to construct \(P\): a
required function Gloc
(of type KINBBDLocalFn
) which approximates the
nonlinear system function \(G(u) \approx F(u)\) and which is computed
locally, and an optional function Gcomm
(of type KINBBDCommFn
) which
performs all interprocess communication necessary to evaluate the approximate
function \(G\). These are in addition to the usersupplied nonlinear system
function that evaluates \(F(u)\). Both functions take as input the same
pointer user_data
as that passed by the user to KINSetUserData()
and
passed to the user’s function func
, and neither function has a return value.
The user is responsible for providing space (presumably within user_data
)
for components of u
that are communicated by Gcomm
from the other
processes, and that are then used by Gloc
, which should not do any
communication.

typedef int (*KINBBDLocalFn)(sunindextype Nlocal, N_Vector u, N_Vector gval, void *user_data)
This
Gloc
function computes \(G(u)\), and outputs the resulting vector asgval
. Arguments:
Nlocal
– is the local vector length.u
– is the current value of the iterate.gval
– is the output vector.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.
 Return value:
An
KINBBDLocalFn
function type should return 0 to indicate success, or nonzero if an error occured. 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 (*KINBBDCommFn)(sunindextype Nlocal, N_Vector u, void *user_data)
This
Gcomm
function performs all interprocessor communications necessary for the execution of theGloc
function above, using the input vectorsu
. Arguments:
Nlocal
– is the local vector length.u
– is the current value of the iterate.user_data
– is a pointer to user data, the same as theuser_data
parameter passed toKINSetUserData()
.
 Return value:
An
KINBBDLocalFn
function type should return 0 to indicate success, or nonzero if an error occured. 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 functionfunc
with the sameu
argument. ThusGcomm
can omit any communications done byfunc
if relevant to the evaluation ofGloc
. If all necessary communication was done infunc
, thenGcomm = NULL
can be passed in the call toKINBBDPrecInit()
.
Besides the header files required for the integration of the DAE problem (see
§8.4.3), to use the KINBBDPRE module, the main program
must include the header file kin_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 §8.4.4 are not bold.
Initialize parallel or multithreaded environment (if appropriate)
Create the SUNDIALS context object
Set the problem dimensions etc.
Create the vector with the initial guess
Create matrix object (if appropriate)
Create linear solver object (if appropriate)
When creating the iterative linear solver object, specify the use of right preconditioning (
SUN_PREC_RIGHT
) as KINSOL only supports right preconditioning.Create nonlinear solver object (if appropriate)
Create KINSOL object
Initialize KINSOL solver
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
KINSetPreconditioner()
optional input function.Initialize the KINBBDPRE preconditioner module
Call
KINBBDPrecInit()
to allocate memory and initialize the internal preconditioner data. The last two arguments ofKINBBDPrecInit()
are the two usersupplied functions described above.Set optional inputs
Solve problem
Get optional outputs
Additional optional outputs associated with KINBBDPRE are available by way of two routines described below,
KINBBDPrecGetWorkSpace()
andKINBBDPrecGetNumGfnEvals()
.Deallocate memory
Finalize MPI, if used
The usercallable functions that initialize or reinitialize the KINBBDPRE preconditioner module are described next.

int KINBBDPrecInit(void *kin_mem, sunindextype Nlocal, sunindextype mudq, sunindexype mldq, sunindextype mukeep, sunindextype mlkeep, realtype dq_rel_u, KINBBDLocalFn Gloc, KINBBDCommFn Gcomm)
The function
KINBBDPrecInit()
initializes and allocates memory for the KINBBDPRE preconditioner. Arguments:
kin_mem
– pointer to the KINSOL memory block.Nlocal
– local vector length.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_u
– the relative increment in components ofu
used in the difference quotient approximations. The default is \(\texttt{dq\_rel\_u} = \sqrt{\text{unit roundoff}}\) , which can be specified by passingdq_rel_u= 0.0
.Gloc
– the CC function which computes the approximation \(G(u) \approx F(u)\).Gcomm
– the optional CC function which performs all interprocess communication required for the computation of \(G(u)\).
 Return value:
KINLS_SUCCESS
– The call toKINBBDPrecInit()
was successful.KINLS_MEM_NULL
– Thekin_mem
pointer wasNULL
.KINLS_MEM_FAIL
– A memory allocation request has failed.KINLS_LMEM_NULL
– The KINLS linear solver interface has not been initialized.KINLS_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 with 0 orNlocal1
accordingly.The halfbandwidths
mudq
andmldq
need not be the true halfbandwidths of the Jacobian of the local block of \(G\), when smaller values may provide greater efficiency.Also, the halfbandwidths
mukeep
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 for every process.
The following two optional output functions are available for use with the KINBBDPRE module:

int KINBBDPrecGetWorkSpace(void *kin_mem, long int *lenrwBBDP, long int *leniwBBDP)
The function
KINBBDPrecGetWorkSpace()
returns the local sizes of the KINBBDPRE real and integer workspaces. Arguments:
kin_mem
– pointer to the KINSOL solver object.lenrwBBDP
– local number of real values in the KINBBDPRE workspace.leniwBBDP
– local number of integer values in the KINBBDPRE workspace.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer wasNULL
.KINLS_PMEM_NULL
– The KINBBDPRE preconditioner has not been initialized.
 Notes:
The workspace requirements reported by this routine correspond only to memory allocated within the KINBBDPRE 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
KINGetLinWorkSpace()
function.

int KINBBDPrecGetNumGfnEvals(void *kin_mem, long int *ngevalsBBDP)
The function
KINBBDPrecGetNumGfnEvals()
returns the cumulative number of calls to the userGres
function due to the finite difference approximation of the Jacobian blocks used within KINBBDPRE’s preconditioner setup function. Arguments:
kin_mem
– pointer to the KINSOL solver object.ngevalsBBDP
– the cumulative number of calls to the userGres
function.
 Return value:
KINLS_SUCCESS
– The optional output value has been successfully set.KINLS_MEM_NULL
– Thekin_mem
pointer wasNULL
.KINLS_PMEM_NULL
– The KINBBDPRE preconditioner has not been initialized.
In addition to the ngevalsBBDP
evaluations of Gres
, the costs associated
with KINBBDPRE also includes nlinsetups
LU factorizations, nlinsetups
calls to Gcomm
, npsolves
banded backsolve calls, and nrevalsLS
residual function evaluations, where nlinsetups
is an optional KINSOL output
(see §8.4.5.5.1), and npsolves
and
nrevalsLS
are linear solver optional outputs (see
§8.4.5.5.2).
8.4.8. Alternative to KINSOL for difficult systems
A nonlinear system \(F(u) = 0\) may be difficult to solve with KINSOL (or any other nonlinear system solver) for a variety of reasons. The possible reasons include high nonlinearity, small region of convergence, and lack of a good initial guess. For systems with such difficulties, there is an alternative approach that may be more successful. This is an old idea, but deserves some new attention.
If the nonlinear system is \(F(u) = 0\), consider instead the ODE system \(du/dt =  M^{1} F(u)\), where \(M\) is a nonsingular matrix that is an approximation (even a crude approximation) to the system Jacobian \(F_u = dF/du\). Whatever \(M\) is, if this ODE is solved until it reaches a steady state \(u^*\), then \(u^*\) is a zero of the righthand side of the ODE, and hence a solution to \(F(u) = 0\). There is no issue of having a close enough initial guess.
A further basis for this choice of ODE is the following: If \(M\) approximates \(F_u\), then the Jacobian of the ODE system, \(M^{1}F_u\), is approximately equal to \(I\) where \(I\) is the identity matrix. This means that (in a local approximation sense) the solution modes of the ODE behave like \(\exp(t)\), and that asymptotically the approach to the steady state goes as \(\exp(t)\). Of course, the closer \(M\) is to \(F_u\), the better this basis applies.
Using (say) CVODE to solve the above ODE system requires, in addition to the objective function \(F(u)\), the calculation of a suitable matrix \(M\) and its inverse, or at least a routine that solves linear systems \(Mx = b\). This is similar to the KINSOL requirement of supplying the system Jacobian \(J\) (or solutions to \(Jx = b\)), but differs in that \(M\) may be simpler than \(J\) and hence easier to deal with. Depending on the nature of \(M\), this may be handled best with a direct solver, or with a preconditioned Krylov solver. The latter calls for the use of a preconditioner \(P\) that may be a crude approximation to \(M\), hence even easier to solve. Note if using ARKODE, the ODE system may be posed in the linearly implicit from \(M du/dt = F(u)\) where \(M\) is the “mass matrix” for the system. This use case requires supplying ARKODE with a function to evaluate \(M\) or to compute its action on a vector (\(Mv = w\)) and attaching a linear solver (direct or iterative) to solve the linear systems \(Mx = b\).
The solution of the ODE may be made easier by solving instead the equivalent DAE, \(M du/dt + F(u) = 0\). Applying IDA to this system requires solutions to linear systems whose matrix is the DAE system Jacobian, \(J = F_u + \alpha M\), where \(\alpha\) is the scalar coefficient \(c_j\) supplied to the user’s Jacobian or preconditioner routine. Selecting a preconditioned Krylov method requires an approximation to this Jacobian as preconditioner \(P\). Given that \(M\) approximates \(F_u\) (possibly crudely), the appropriate approximation to \(J\) is \(P = M + \alpha M = (1 + \alpha)M\). Again the user must supply a routine that solves linear systems \(Px = b\), or \(M x = b/(1 + \alpha)\). If M is too difficult to solve, than an approximation \(M'\) that is easier can be substituted, as long as it achieves convergence. As always, there is a tradeoff between the expense of solving \(M\) and the difficulty of achieving convergence in the linear solver.
For the solution of either the ODE or DAE system above, the chances for convergence can be improved with a piecewise constant choice for \(M\). Specifically, starting from an initial guess \(u_0\), an initial choice for \(M\) might be \(M_0 = F_u(u_0)\), or some approximation to this Jacobian. Then one could integrate \(M_0 du/dt + F(u) = 0\) from \(t = 0\) to \(t = T\) for some sizable value \(T\), evaluate \(F_u(u(T))\), and take \(M_1\) to be an approximation to that Jacobian. Then integrate using \(M_1\) from \(t = T\) to \(t = 2T\), and repeat the process until it converges to a steady state.