5.4.4. Using CVODES for Adjoint Sensitivity Analysis

This chapter describes the use of CVODES to compute sensitivities of derived functions using adjoint sensitivity analysis. As mentioned before, the adjoint sensitivity module of CVODES provides the infrastructure for integrating backward in time any system of ODEs that depends on the solution of the original IVP, by providing various interfaces to the main CVODES integrator, as well as several supporting user-callable functions. For this reason, in the following sections we refer to the backward problem and not to the adjoint problem when discussing details relevant to the ODEs that are integrated backward in time. The backward problem can be the adjoint problem (5.21) or (5.24), and can be augmented with some quadrature differential equations.

CVODES uses various constants for both input and output. These are defined as needed in this chapter, but for convenience are also listed separately in §5.5.

We begin with a brief overview, in the form of a skeleton user program. Following that are detailed descriptions of the interface to the various user-callable functions and of the user-supplied functions that were not already described in §5.4.1.

5.4.4.1. A skeleton of the user’s main program

The following is a skeleton of the user’s main program as an application of CVODES. The user program is to have these steps in the order indicated, unless otherwise noted. For the sake of brevity, we defer many of the details to the later sections. As in §5.4.1.2, most steps are independent of the N_Vector, SUNMatrix, SUNLinearSolver, and SUNNonlinearSolver implementations used. For the steps that are not, refer to Chapters §9, §10, §11, and §12 for the specific name of the function to be called or macro to be referenced.

Steps that are unchanged from the skeleton programs presented in §5.4.1.2, §5.4.3.1, and §5.4.2 are grayed out and new or modified steps are in bold.

  1. Initialize parallel or multi-threaded environment, if appropriate

  2. Create the SUNDIALS context object

  3. Set initial conditions for the forward problem

  4. Create CVODES object for the forward problem

  5. Initialize CVODES for the forward problem

  6. Specify integration tolerances for forward problem

  7. Create matrix object for the forward problem

  8. Create linear solver object for the forward problem

  9. Set linear solver optional inputs for the forward problem

  10. Attach linear solver module for the forward problem

  11. Set optional inputs for the forward problem

  12. Create nonlinear solver object for the forward problem

  13. Attach nonlinear solver module for the forward problem

  14. Set nonlinear solver optional inputs for the forward problem

  15. Initialize quadrature problem or problems for forward problems

  16. Initialize forward sensitivity problem

  17. Specify rootfinding

  18. Allocate space for the adjoint computation

    Call CVodeAdjInit() to allocate memory for the combined forward-backward problem. This call requires Nd, the number of steps between two consecutive checkpoints. CVodeAdjInit() also specifies the type of interpolation used (see §5.2.9).

  19. Integrate forward problem

    Call CVodeF(), a wrapper for the CVODES main integration function CVode(), either in CV_NORMAL mode to the time tout or in CV_ONE_STEP mode inside a loop (if intermediate solutions of the forward problem are desired). The final value of tret is then the maximum allowable value for the endpoint \(T\) of the backward problem.

  20. Set problem dimensions etc. for the backward problem

    This generally includes the backward problem vector length NB, and possibly the local vector length NBlocal.

  21. Set initial values for the backward problem

    Set the endpoint time tB0 = T, and set the corresponding vector yB0 at which the backward problem starts.

  22. Create the backward problem

    Call CVodeCreateB(), a wrapper for CVodeCreate(), to create the CVODES memory block for the new backward problem. Unlike CVodeCreate(), the function CVodeCreateB() does not return a pointer to the newly created memory block. Instead, this pointer is attached to the internal adjoint memory block (created by CVodeAdjInit()) and returns an identifier called which that the user must later specify in any actions on the newly created backward problem.

  23. Allocate memory for the backward problem

    Call CVodeInitB() (or CVodeInitBS(), when the backward problem depends on the forward sensitivities). The two functions are actually wrappers for CVodeInit() and allocate internal memory, specify problem data, and initialize CVODES at tB0 for the backward problem.

  24. Specify integration tolerances for backward problem

    Call CVodeSStolerancesB() or CVodeSVtolerancesB() to specify a scalar relative tolerance and scalar absolute tolerance or scalar relative tolerance and a vector of absolute tolerances, respectively. The functions are wrappers for CVodeSStolerances() and CVodeSVtolerances(), but they require an extra argument which, the identifier of the backward problem returned by CVodeCreateB().

  25. Create matrix object for the backward problem

    If a nonlinear solver requiring a linear solve will be used (e.g., the the default Newton iteration) and the linear solver will be a direct linear solver, then a template Jacobian matrix must be created by calling the appropriate constructor function defined by the particular SUNMatrix implementation.

    For the native SUNDIALS SUNMatrix implementations, the matrix object may be created using a call of the form SUN***Matrix(...) where *** is the name of the matrix (see §10 for details).

  26. Create linear solver object for the backward problem

    If a nonlinear solver requiring a linear solver is chosen (e.g., the default Newton iteration), then the desired linear solver object for the backward problem must be created by calling the appropriate constructor function defined by the particular SUNLinearSolver implementation.

    For any of the SUNDIALS-supplied SUNLinearSolver implementations, the linear solver object may be created using a call of the form

    SUNLinearSolver LS = SUNLinSol_*(...);

    where * can be replaced with “Dense”, “SPGMR”, or other options, as discussed in §5.4.1.3.5 and Chapter §11.

    Note that it is not required to use the same linear solver module for both the forward and the backward problems; for example, the forward problem could be solved with the SUNLINSOL_BAND linear solver module and the backward problem with SUNLINSOL_SPGMR linear solver module.

  27. Set linear solver interface optional inputs for the backward problem

    Call *Set* functions from the selected linear solver module to change optional inputs specific to that linear solver. See the documentation for each SUNLinearSolver module in Chapter §11.

  28. Attach linear solver module for the backward problem

    If a nonlinear solver requiring a linear solver is chosen for the backward problem (e.g., the default Newton iteration), then initialize the CVLS linear solver interface by attaching the linear solver object (and matrix object, if applicable) with the call to CVodeSetLinearSolverB()

    Alternately, if the CVODES-specific diagonal linear solver module, CVDIAG, is desired, initialize the linear solver module and attach it to CVODES with a call to CVDiagB().

  29. Set optional inputs for the backward problem

    Call CVodeSet*B functions to change from their default values any optional inputs that control the behavior of CVODES. Unlike their counterparts for the forward problem, these functions take an extra argument which, the identifier of the backward problem returned by CVodeCreateB().

  30. Create nonlinear solver object for the backward problem (optional)

    If using a non-default nonlinear solver for the backward problem, then create the desired nonlinear solver object by calling the appropriate constructor function defined by the particular SUNNonlinearSolver implementation (e.g., NLSB = SUNNonlinSol_***(...); where *** is the name of the nonlinear solver.

  31. Attach nonlinear solver module for the backward problem (optional)

    If using a non-default nonlinear solver for the backward problem, then initialize the nonlinear solver interface by attaching the nonlinear solver object by calling CVodeSetNonlinearSolverB().

  32. Initialize quadrature calculation

    If additional quadrature equations must be evaluated, call CVodeQuadInitB() or CVodeQuadInitBS() (if quadrature depends also on the forward sensitivities). These functions are wrappers around CVodeQuadInit() and can be used to initialize and allocate memory for quadrature integration. Optionally, call CVodeSetQuad*B functions to change from their default values optional inputs that control the integration of quadratures during the backward phase.

  33. Integrate backward problem

    Call CVodeB(), a second wrapper around the CVODES main integration function CVode(), to integrate the backward problem from tB0. This function can be called either in CV_NORMAL or CV_ONE_STEP mode. Typically, CVodeB() will be called in CV_NORMAL mode with an end time equal to the initial time \(t_0\) of the forward problem.

  34. Extract quadrature variables

    If applicable, call CVodeGetQuadB(), a wrapper around CVodeGetQuad(), to extract the values of the quadrature variables at the time returned by the last call to CVodeB().

  35. Deallocate memory

    Upon completion of the backward integration, call all necessary deallocation functions. These include appropriate destructors for the vectors y and yB, a call to CVodeFree() to free the CVODES memory block for the forward problem. If one or more additional Adjoint Sensitivity Analyses are to be done for this problem, a call to CVodeAdjFree() may be made to free and deallocate memory allocated for the backward problems, followed by a call to CVodeAdjInit().

    Free the nonlinear solver memory for the forward and backward problems

    Free linear solver and matrix memory for the forward and backward problems

  36. Finalize MPI, if used

The above user interface to the adjoint sensitivity module in CVODES was motivated by the desire to keep it as close as possible in look and feel to the one for ODE IVP integration. Note that if steps back_start-back_end are not present, a program with the above structure will have the same functionality as one described in §5.4.1.2 for integration of ODEs, albeit with some overhead due to the checkpointing scheme.

If there are multiple backward problems associated with the same forward problem, repeat steps back_start-back_end above for each successive backward problem. In the process, each call to CVodeCreateB() creates a new value of the identifier which.

5.4.4.2. User-callable functions for adjoint sensitivity analysis

5.4.4.2.1. Adjoint sensitivity allocation and deallocation functions

After the setup phase for the forward problem, but before the call to CVodeF(), memory for the combined forward-backward problem must be allocated by a call to the function CVodeAdjInit(). The form of the call to this function is

int CVodeAdjInit(void *cvode_mem, long int Nd, int interpType)

The function CVodeAdjInit() updates CVODES memory block by allocating the internal memory needed for backward integration. Space is allocated for the Nd = N_d interpolation data points, and a linked list of checkpoints is initialized.

Arguments:
  • cvode_mem – is the pointer to the CVODES memory block returned by a previous call to CVodeCreate().

  • Nd – is the number of integration steps between two consecutive checkpoints.

  • interpType – specifies the type of interpolation used and can be CV_POLYNOMIAL or CV_HERMITE , indicating variable-degree polynomial and cubic Hermite interpolation, respectively see §5.2.9.

Return value:
  • CV_SUCCESSCVodeAdjInit() was successful.

  • CV_MEM_FAIL – A memory allocation request has failed.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_ILL_INPUT – One of the parameters was invalid: Nd was not positive or interpType is not one of the CV_POLYNOMIAL or CV_HERMITE.

Notes:

The user must set Nd so that all data needed for interpolation of the forward problem solution between two checkpoints fits in memory. CVodeAdjInit() attempts to allocate space for 2*Nd+3 variables of type N_Vector. If an error occurred, CVodeAdjInit() also sends a message to the error handler function.

int CVodeAdjReInit(void *cvode_mem)

The function CVodeAdjReInit() reinitializes the CVODES memory block for ASA, assuming that the number of steps between check points and the type of interpolation remain unchanged.

Arguments:
  • cvode_mem – is the pointer to the CVODES memory block returned by a previous call to CVodeCreate().

Return value:
  • CV_SUCCESSCVodeAdjReInit() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() was not previously called.

Notes:

The list of check points (and associated memory) is deleted. The list of backward problems is kept. However, new backward problems can be added to this list by calling CVodeCreateB(). If a new list of backward problems is also needed, then free the adjoint memory (by calling CVodeAdjFree()) and reinitialize ASA with CVodeAdjInit(). The CVODES memory for the forward and backward problems can be reinitialized separately by calling CVodeReInit() and CVodeReInitB(), respectively.

void CVodeAdjFree(void *cvode_mem)

The function CVodeAdjFree() frees the memory related to backward integration allocated by a previous call to CVodeAdjInit().

Argument:
  • cvode_mem – is the pointer to the CVODES memory block returned by a previous call to CVodeCreate().

Return value:

The function has no return value.

Notes:

This function frees all memory allocated by CVodeAdjInit(). This includes workspace memory, the linked list of checkpoints, memory for the interpolation data, as well as the CVODES memory for the backward integration phase. Unless one or more further calls to CVodeAdjInit() are to be made, CVodeAdjFree() should not be called by the user, as it is invoked automatically by CVodeFree().

5.4.4.2.2. Forward integration function

The function CVodeF() is very similar to the CVODES function CVode() in that it integrates the solution of the forward problem and returns the solution in y. At the same time, however, CVodeF() stores checkpoint data every Nd integration steps. CVodeF() can be called repeatedly by the user. Note that CVodeF() is used only for the forward integration pass within an Adjoint Sensitivity Analysis. It is not for use in Forward Sensitivity Analysis; for that, see §5.4.3. The call to this function has the form

int CVodeF(void *cvode_mem, sunrealtype tout, N_Vector yret, sunrealtype *tret, int itask, int *ncheck)

The function CVodeF() integrates the forward problem over an interval in \(t\) and saves checkpointing data.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • tout – the next time at which a computed solution is desired.

  • yret – the computed solution vector \(y\).

  • tret – the time reached by the solver output.

  • itask – output mode a flag indicating the job of the solver for the next step. The CV_NORMAL task is to have the solver take internal steps until it has reached or just passed the user-specified tout parameter. The solver then interpolates in order to return an approximate value of \(y(\text{tout})\). The CV_ONE_STEP option tells the solver to just take one internal step and return the solution at the point reached by that step.

  • ncheck – the number of internal checkpoints stored so far.

Return value:
  • CV_SUCCESSCVodeF() succeeded.

  • CV_TSTOP_RETURNCVodeF() succeeded by reaching the optional stopping point.

  • CV_ROOT_RETURNCVodeF() succeeded and found one or more roots. In this case, tret is the location of the root. If nrtfn > 1 , call CVodeGetRootInfo() to see which \(g_i\) were found to have a root.

  • CV_NO_MALLOC – The function CVodeInit() has not been previously called.

  • CV_ILL_INPUT – One of the inputs to CVodeF() is illegal.

  • CV_TOO_MUCH_WORK – The solver took mxstep internal steps but could not reach tout.

  • CV_TOO_MUCH_ACC – The solver could not satisfy the accuracy demanded by the user for some internal step.

  • CV_ERR_FAILURE – Error test failures occurred too many times during one internal time step or occurred with \(|h| = h_{min}\).

  • CV_CONV_FAILURE – Convergence test failures occurred too many times during one internal time step or occurred with \(|h| = h_{min}\).

  • CV_LSETUP_FAIL – The linear solver’s setup function failed in an unrecoverable manner.

  • CV_LSOLVE_FAIL – The linear solver’s solve function failed in an unrecoverable manner.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_MEM_FAIL – A memory allocation request has failed in an attempt to allocate space for a new checkpoint.

Notes:

All failure return values are negative and therefore a test flag\(< 0\) will trap all CVodeF() failures. At this time, CVodeF() stores checkpoint information in memory only. Future versions will provide for a safeguard option of dumping checkpoint data into a temporary file as needed. The data stored at each checkpoint is basically a snapshot of the CVODES internal memory block and contains enough information to restart the integration from that time and to proceed with the same step size and method order sequence as during the forward integration. In addition, CVodeF() also stores interpolation data between consecutive checkpoints so that, at the end of this first forward integration phase, interpolation information is already available from the last checkpoint forward. In particular, if no checkpoints were necessary, there is no need for the second forward integration phase.

Warning

It is illegal to change the integration tolerances between consecutive calls to CVodeF(), as this information is not captured in the checkpoint data.

5.4.4.2.3. Backward problem initialization functions

The functions CVodeCreateB() and CVodeInitB() (or CVodeInitBS()) must be called in the order listed. They instantiate a CVODES solver object, provide problem and solution specifications, and allocate internal memory for the backward problem.

int CVodeCreateB(void *cvode_mem, int lmmB, int *which)

The function CVodeCreateB() instantiates a CVODES solver object and specifies the solution method for the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • lmmB – specifies the linear multistep method and may be one of two possible values: CV_ADAMS or CV_BDF.

  • which – contains the identifier assigned by CVODES for the newly created backward problem. Any call to CVode*B functions requires such an identifier.

Return value:
  • CV_SUCCESS – The call to CVodeCreateB() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_MEM_FAIL – A memory allocation request has failed.

There are two initialization functions for the backward problem – one for the case when the backward problem does not depend on the forward sensitivities, and one for the case when it does. These two functions are described next.

int CVodeInitB(void *cvode_mem, int which, CVRhsFnB rhsB, sunrealtype tB0, N_Vector yB0)

The function CVodeInitB() provides problem specification, allocates internal memory, and initializes the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • rhsB – is the CVRhsFnB function which computes \(f_B\) , the right-hand side of the backward ODE problem.

  • tB0 – specifies the endpoint \(T\) where final conditions are provided for the backward problem, normally equal to the endpoint of the forward integration.

  • yB0 – is the initial value at \(t =\) tB0 of the backward solution.

Return value:
  • CV_SUCCESS – The call to CVodeInitB() was successful.

  • CV_NO_MALLOC – The function CVodeInit() has not been previously called.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_BAD_TB0 – The final time tB0 was outside the interval over which the forward problem was solved.

  • CV_ILL_INPUT – The parameter which represented an invalid identifier, or either yB0 or rhsB was NULL.

Notes:

The memory allocated by CVodeInitB() is deallocated by the function CVodeAdjFree().

The function CVodeInitB() initializes the backward problem when it does not depend on the forward sensitivities. It is essentially a wrapper for CVodeInit() with some particularization for backward integration, as described below.

For the case when backward problem also depends on the forward sensitivities, user must call CVodeInitBS() instead of CVodeInitB(). Only the third argument of each function differs between these two functions.

int CVodeInitBS(void *cvode_mem, int which, CVRhsFnBS rhsBS, sunrealtype tB0, N_Vector yB0)

The function CVodeInitBS() provides problem specification, allocates internal memory, and initializes the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • rhsBS – is the CVRhsFnBS function which computes \(f_B\) , the right-hand side of the backward ODE problem.

  • tB0 – specifies the endpoint \(T\) where final conditions are provided for the backward problem.

  • yB0 – is the initial value at \(t =\) tB0 of the backward solution.

Return value:
  • CV_SUCCESS – The call to CVodeInitB() was successful.

  • CV_NO_MALLOC – The function CVodeInit() has not been previously called.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_BAD_TB0 – The final time tB0 was outside the interval over which the forward problem was solved.

  • CV_ILL_INPUT – The parameter which represented an invalid identifier, either yB0 or rhsBS was NULL , or sensitivities were not active during the forward integration.

Notes:

The memory allocated by CVodeInitBS() is deallocated by the function CVodeAdjFree().

The function CVodeReInitB() reinitializes CVODES for the solution of a series of backward problems, each identified by a value of the parameter which. CVodeReInitB() is essentially a wrapper for CVodeReInit(), and so all details given for CVodeReInit() apply here. Also note that CVodeReInitB() can be called to reinitialize the backward problem even it has been initialized with the sensitivity-dependent version CVodeInitBS(). Before calling CVodeReInitB() for a new backward problem, call any desired solution extraction functions CVodeGet** associated with the previous backward problem. The call to the CVodeReInitB() function has the form

int CVodeReInitB(void *cvode_mem, int which, sunrealtype tB0, N_Vector yB0)

The function CVodeReInitB() reinitializes a CVODES backward problem.

Arguments:
  • cvode_mem – pointer to CVODES memory block returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • tB0 – specifies the endpoint \(T\) where final conditions are provided for the backward problem.

  • yB0 – is the initial value at \(t =\) tB0 of the backward solution.

Return value:
  • CV_SUCCESS – The call to CVodeReInitB() was successful.

  • CV_NO_MALLOC – The function CVodeInit() has not been previously called.

  • CV_MEM_NULL – The cvode_mem memory block pointer was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_BAD_TB0 – The final time tB0 is outside the interval over which the forward problem was solved.

  • CV_ILL_INPUT – The parameter which represented an invalid identifier, or yB0 was NULL.

5.4.4.2.4. Tolerance specification functions for backward problem

One of the following two functions must be called to specify the integration tolerances for the backward problem. Note that this call must be made after the call to CVodeInitB() or CVodeInitBS().

int CVodeSStolerancesB(void *cvode_mem, int which, sunrealtype reltolB, sunrealtype abstolB)

The function CVodeSStolerancesB() specifies scalar relative and absolute tolerances.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • reltolB – is the scalar relative error tolerance.

  • abstolB – is the scalar absolute error tolerance.

Return value:
  • CV_SUCCESS – The call to CVodeSStolerancesB() was successful.

  • CV_MEM_NULL – The CVODES memory block was not initialized through a previous call to CVodeCreate().

  • CV_NO_MALLOC – The allocation function CVodeInit() has not been called.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_ILL_INPUT – One of the input tolerances was negative.

int CVodeSVtolerancesB(void *cvode_mem, int which, reltolBabstolB)

The function CVodeSVtolerancesB() specifies scalar relative tolerance and vector absolute tolerances.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • reltol – is the scalar relative error tolerance.

  • abstol – is the vector of absolute error tolerances.

Return value:
  • CV_SUCCESS – The call to CVodeSVtolerancesB() was successful.

  • CV_MEM_NULL – The CVODES memory block was not initialized through a previous call to CVodeCreate().

  • CV_NO_MALLOC – The allocation function CVodeInit() has not been called.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_ILL_INPUT – The relative error tolerance was negative or the absolute tolerance had a negative component.

Notes:

This choice of tolerances is important when the absolute error tolerance needs to be different for each component of the state vector \(y\).

5.4.4.2.5. Linear solver initialization functions for backward problem

All CVODES linear solver modules available for forward problems are available for the backward problem. They should be created as for the forward problem and then attached to the memory structure for the backward problem using the following functions.

int CVodeSetLinearSolverB(void *cvode_mem, int which, SUNLinearSolver LS, SUNMatrix A)

The function CVodeSetLinearSolverB() attaches a generic SUNLinearSolver object LS and corresponding template Jacobian SUNMatrix object A to CVODES, initializing the CVLS linear solver interface for solution of the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – represents the identifier of the backward problem returned by CVodeCreateB().

  • LS – SUNLINSOL object to use for solving linear systems for the backward problem.

  • A – SUNMATRIX object for used as a template for the Jacobian for the backward problem or NULL if not applicable.

Return value:
  • CVLS_SUCCESS – The CVLS initialization was successful.

  • CVLS_MEM_NULL – The cvode_mem pointer is NULL.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

  • CVLS_MEM_FAIL – A memory allocation request failed.

  • CVLS_NO_ADJ – The function CVAdjInit has not been previously called.

Notes:

If LS is a matrix-based linear solver, then the template Jacobian matrix J will be used in the solve process, so if additional storage is required within the SUNMatrix object (e.g., for factorization of a banded matrix), ensure that the input object is allocated with sufficient size (see the documentation of the particular SUNMatrix type in §10). The previous routines CVDlsSetLinearSolverB and CVSpilsSetLinearSolverB are now wrappers for this routine, and may still be used for backward-compatibility. However, these will be deprecated in future releases, so we recommend that users transition to the new routine name soon.

int CVDiagB(void *cvode_mem, int which)

The function CVDiagB selects the CVDIAG linear solver for the solution of the backward problem. The user’s main program must include the cvodes_diag.h header file.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – represents the identifier of the backward problem returned by CVodeCreateB().

Return value:
  • CVDIAG_SUCCESS – The CVDIAG initialization was successful.

  • CVDIAG_MEM_NULL – The cvode_mem pointer is NULL.

  • CVDIAG_ILL_INPUT – The CVDIAG solver is not compatible with the current NVECTOR module.

  • CVDIAG_MEM_FAIL – A memory allocation request failed.

Notes:

The CVDIAG solver is the simplest of all of the available CVODES linear solver interfaces. The CVDIAG solver uses an approximate diagonal Jacobian formed by way of a difference quotient. The user does not have the option of supplying a function to compute an approximate diagonal Jacobian.

5.4.4.2.6. Nonlinear solver initialization function for backward problem

All CVODES nonlinear solver modules available for forward problems are available for the backward problem. As with the forward problem CVODES uses the SUNNonlinearSolver implementation of Newton’s method defined by the SUNNONLINSOL_NEWTON module by default.

To specify a different nonlinear solver for the backward problem, the user’s program must create a SUNNonlinearSolver object by calling the appropriate constructor routine. The user must then attach the SUNNonlinearSolver object by calling CVodeSetNonlinearSolverB(), as documented below.

When changing the nonlinear solver in CVODES, CVodeSetNonlinearSolverB() must be called after CVodeInitB(). If any calls to CVodeB() have been made, then CVODES will need to be reinitialized by calling CVodeReInitB() to ensure that the nonlinear solver is initialized correctly before any subsequent calls to CVodeB().

int CVodeSetNonlinearSolverB(void *cvode_mem, int which, SUNNonlinearSolver NLS)

The function CVodeSetNonLinearSolverB() attaches a SUNNONLINEARSOLVER object (NLS) to CVODES for the solution of the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – represents the identifier of the backward problem returned by CVodeCreateB().

  • NLS – SUNNONLINSOL object to use for solving nonlinear systems for the backward problem.

Return value:
  • CV_SUCCESS – The nonlinear solver was successfully attached.

  • CV_MEM_NULL – The cvode_mem pointer is NULL.

  • CVLS_NO_ADJ – The function CVAdjInit has not been previously called.

  • CV_ILL_INPUT – The parameter which represented an invalid identifier or the SUNNONLINSOL object is NULL , does not implement the required nonlinear solver operations, is not of the correct type, or the residual function, convergence test function, or maximum number of nonlinear iterations could not be set.

5.4.4.2.7. Backward integration function

The function CVodeB() performs the integration of the backward problem. It is essentially a wrapper for the CVODES main integration function CVode() and, in the case in which checkpoints were needed, it evolves the solution of the backward problem through a sequence of forward-backward integration pairs between consecutive checkpoints. The first run of each pair integrates the original IVP forward in time and stores interpolation data; the second run integrates the backward problem backward in time and performs the required interpolation to provide the solution of the IVP to the backward problem.

The function CVodeB() does not return the solution yB itself. To obtain that, call the function CVodeGetB(), which is also described below.

The CVodeB() function does not support rootfinding, unlike CVodeF(), which supports the finding of roots of functions of \((t,y)\). If rootfinding was performed by CVodeF(), then for the sake of efficiency, it should be disabled for CVodeB() by first calling CVodeRootInit() with nrtfn = 0.

The call to CVodeB() has the form

int CVodeB(void *cvode_mem, sunrealtype tBout, int itaskB)

The function CVodeB() integrates the backward ODE problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • tBout – the next time at which a computed solution is desired.

  • itaskB – output mode a flag indicating the job of the solver for the next step. The CV_NORMAL task is to have the solver take internal steps until it has reached or just passed the user-specified value tBout. The solver then interpolates in order to return an approximate value of \(yB(\texttt{tBout})\). The CV_ONE_STEP option tells the solver to take just one internal step in the direction of tBout and return.

Return value:
  • CV_SUCCESSCVodeB() succeeded.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_NO_BCK – No backward problem has been added to the list of backward problems by a call to CVodeCreateB().

  • CV_NO_FWD – The function CVodeF() has not been previously called.

  • CV_ILL_INPUT – One of the inputs to CVodeB() is illegal.

  • CV_BAD_ITASK – The itaskB argument has an illegal value.

  • CV_TOO_MUCH_WORK – The solver took mxstep internal steps but could not reach tBout.

  • CV_TOO_MUCH_ACC – The solver could not satisfy the accuracy demanded by the user for some internal step.

  • CV_ERR_FAILURE – Error test failures occurred too many times during one internal time step.

  • CV_CONV_FAILURE – Convergence test failures occurred too many times during one internal time step.

  • CV_LSETUP_FAIL – The linear solver’s setup function failed in an unrecoverable manner.

  • CV_SOLVE_FAIL – The linear solver’s solve function failed in an unrecoverable manner.

  • CV_BCKMEM_NULL – The solver memory for the backward problem was not created with a call to CVodeCreateB().

  • CV_BAD_TBOUT – The desired output time tBout is outside the interval over which the forward problem was solved.

  • CV_REIFWD_FAIL – Reinitialization of the forward problem failed at the first checkpoint corresponding to the initial time of the forward problem.

  • CV_FWD_FAIL – An error occurred during the integration of the forward problem.

Notes:

All failure return values are negative and therefore a test flag < 0 will trap all CVodeB() failures. In the case of multiple checkpoints and multiple backward problems, a given call to CVodeB() in CV_ONE_STEP mode may not advance every problem one step, depending on the relative locations of the current times reached. But repeated calls will eventually advance all problems to tBout.

In the case of multiple checkpoints and multiple backward problems, a given call to CVodeB() in CV_ONE_STEP mode may not advance every problem one step, depending on the relative locations of the current times reached. But repeated calls will eventually advance all problems to tBout.

To obtain the solution yB to the backward problem, call the function CVodeGetB() as follows:

int CVodeGetB(void *cvode_mem, int which, sunrealtype *tret, N_Vector yB)

The function CVodeGetB() provides the solution yB of the backward ODE problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • which – the identifier of the backward problem.

  • tret – the time reached by the solver output.

  • yB – the backward solution at time tret.

Return value:
  • CV_SUCCESSCVodeGetB() was successful.

  • CV_MEM_NULLcvode_mem is NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_ILL_INPUT – The parameter which is an invalid identifier.

Warning

The user must allocate space for yB. To obtain the solution associated with a given backward problem at some other time within the last integration step, first obtain a pointer to the proper CVODES memory structure by calling CVodeGetAdjCVodeBmem() and then use it to call CVodeGetDky().

5.4.4.2.8. Adjoint sensitivity optional input

At any time during the integration of the forward problem, the user can disable the checkpointing of the forward sensitivities by calling the following function:

int CVodeAdjSetNoSensi(void *cvode_mem)

The function CVodeAdjSetNoSensi() instructs CVodeF() not to save checkpointing data for forward sensitivities anymore.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

Return value:
  • CV_SUCCESS – The call to CVodeCreateB() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

5.4.4.2.9. Optional input functions for the backward problem

As for the forward problem there are numerous optional input parameters that control the behavior of the CVODES solver for the backward problem. CVODES provides functions that can be used to change these optional input parameters from their default values which are then described in detail in the remainder of this section, beginning with those for the main CVODES solver and continuing with those for the linear solver interfaces. Note that the diagonal linear solver module has no optional inputs. For the most casual use of CVODES, the reader can skip to §5.4.4.3.

We note that, on an error return, all of the optional input functions send an error message to the error handler function. All error return values are negative, so the test flag < 0 will catch all errors. Finally, a call to a CVodeSet***B function can be made from the user’s calling program at any time and, if successful, takes effect immediately.

5.4.4.2.9.1. Main solver optional input functions

The adjoint module in CVODES provides wrappers for most of the optional input functions defined in §5.4.1.3.10.1. The only difference is that the user must specify the identifier which of the backward problem within the list managed by CVODES.

The optional input functions defined for the backward problem are:

flag = CVodeSetUserDataB(cvode_mem, which, user_dataB);
flag = CVodeSetMaxOrdB(cvode_mem, which, maxordB);
flag = CVodeSetMaxNumStepsB(cvode_mem, which, mxstepsB);
flag = CVodeSetInitStepB(cvode_mem, which, hinB)
flag = CVodeSetMinStepB(cvode_mem, which, hminB);
flag = CVodeSetMaxStepB(cvode_mem, which, hmaxB);
flag = CVodeSetStabLimDetB(cvode_mem, which, stldetB);
flag = CVodeSetConstraintsB(cvode_mem, which, constraintsB);

Their return value flag (of type int) can have any of the return values of their counterparts, but it can also be CV_NO_ADJ if CVodeAdjInit() has not been called, or CV_ILL_INPUT if which was an invalid identifier.

5.4.4.2.9.2. Linear solver interface optional input functions

When using matrix-based linear solver modules, the CVLS solver interface needs a function to compute an approximation to the Jacobian matrix or the linear system for the backward problem. The function to evaluate the Jacobian can be attached through a call to either CVodeSetJacFnB() or CVodeSetJacFnBS(), with the second used when the backward problem depends on the forward sensitivities.

int CVodeSetJacFnB(void *cvode_mem, int which, CVLsJacFnB jacB)

The function CVodeSetJacFnB() specifies the Jacobian approximation function to be used for the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • jacB – user-defined Jacobian approximation function.

Return value:
  • CVLS_SUCCESSCVodeSetJacFnB() succeeded.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_LMEM_NULL – The linear solver has not been initialized with a call to CVodeSetLinearSolverB().

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The previous routine CVDlsSetJacFnB is now deprecated.

int CVodeSetJacFnBS(void *cvode_mem, int which, CVLsJacFnBS jacBS)

The function CVodeSetJacFnBS() specifies the Jacobian approximation function to be used for the backward problem, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • jacBS – user-defined Jacobian approximation function.

Return value:
  • CVLS_SUCCESSCVodeSetJacFnBS() succeeded.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_LMEM_NULL – The linear solver has not been initialized with a call to CVodeSetLinearSolverB().

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The previous routine CVDlsSetJacFnBS is now deprecated.

int CVodeSetLinSysFnB(void *cvode_mem, int which, CVLsLinSysFnB linsysB)

The function CVodeSetLinSysFnB() specifies the linear system approximation function to be used for the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • linsysB – user-defined linear system approximation function.

Return value:
  • CVLS_SUCCESSCVodeSetLinSysFnB() succeeded.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_LMEM_NULL – The linear solver has not been initialized with a call to CVodeSetLinearSolverB().

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

int CVodeSetLinSysFnBS(void *cvode_mem, int which, CVLsLinSysFnBS linsysBS)

The function CVodeSetLinSysFnBS() specifies the linear system approximation function to be used for the backward problem, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • cvode_mem – pointer to the CVODES memory returned by CVodeCreate().

  • which – represents the identifier of the backward problem.

  • linsysBS – user-defined linear system approximation function.

Return value:
  • CVLS_SUCCESSCVodeSetLinSysFnBS() succeeded.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_LMEM_NULL – The linear solver has not been initialized with a call to CVodeSetLinearSolverB().

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

The function CVodeSetLinearSolutionScalingB() can be used to enable or disable solution scaling when using a matrix-based linear solver.

int CVodeSetLinearSolutionScalingB(void *cvode_mem, int which, sunbooleantype onoffB)

The function CVodeSetLinearSolutionScalingB() enables or disables scaling the linear system solution to account for a change in \(\gamma\) in the linear system in the backward problem. For more details see §11.4.1.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – represents the identifier of the backward problem.

  • onoffB – flag to enable SUNTRUE or disable SUNFALSE scaling

Return value:
  • CVLS_SUCCESS – The flag value has been successfully set.

  • CVLS_MEM_NULL – The cvode_mem pointer is NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver interface has not been initialized.

  • CVLS_ILL_INPUT – The attached linear solver is not matrix-based or the linear multistep method type is not BDF.

Notes:

By default scaling is enabled with matrix-based linear solvers when using BDF methods.

int CVodeSetJacTimesB(void *cvode_mem, int which, CVLsJacTimesSetupFnB jsetupB, CVLsJacTimesVecFnB jtvB)

The function CVodeSetJacTimesB() specifies the Jacobian-vector setup and product functions to be used.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • jtsetupB – user-defined function to set up the Jacobian-vector product. Pass NULL if no setup is necessary.

  • jtvB – user-defined Jacobian-vector product function.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The previous routine CVSpilsSetJacTimesB is now deprecated.

int CVodeSetJacTimesBS(void *cvode_mem, int which, CVLsJacTimesVecFnBS jtvBS)

The function CVodeSetJacTimesBS() specifies the Jacobian-vector setup and product functions to be used, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • jtsetupBS – user-defined function to set up the Jacobian-vector product. Pass NULL if no setup is necessary.

  • jtvBS – user-defined Jacobian-vector product function.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The previous routine CVSpilsSetJacTimesBS is now deprecated.

When using the internal difference quotient the user may optionally supply an alternative right-hand side function for use in the Jacobian-vector product approximation for the backward problem by calling CVodeSetJacTimesRhsFnB(). The alternative right-hand side function should compute a suitable (and differentiable) approximation to the right-hand side function provided to CVodeInitB() or CVodeInitBS(). For example, as done in [46] for a forward integration without sensitivity analysis, the alternative function may use lagged values when evaluating a nonlinearity in the right-hand side to avoid differencing a potentially non-differentiable factor.

int CVodeSetJacTimesRhsFnB(void *cvode_mem, int which, CVRhsFn jtimesRhsFn)

The function CVodeSetJacTimesRhsFn() specifies an alternative ODE right-hand side function for use in the internal Jacobian-vector product difference quotient approximation.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • jtimesRhsFn – is the CC function which computes the alternative ODE right-hand side function to use in Jacobian-vector product difference quotient approximations.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULL – The cvode_mem pointer is NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier or the internal difference quotient approximation is disabled.

Notes:

The default is to use the right-hand side function provided to CVodeInit() in the internal difference quotient. If the input right-hand side function is NULL, the default is used. This function must be called after the CVLS linear solver interface has been initialized through a call to CVodeSetLinearSolverB().

int CVodeSetPreconditionerB(void *cvode_mem, int which, CVLPrecSetupFnB psetupB, CVLsPrecSolveFnB psolveB)

The function CVodeSetPrecSolveFnB() specifies the preconditioner setup and solve functions for the backward integration.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • psetupB – user-defined preconditioner setup function.

  • psolveB – user-defined preconditioner solve function.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The psetupB argument may be NULL if no setup operation is involved in the preconditioner. The previous routine CVSpilsSetPrecSolveFnB is now deprecated.

int CVodeSetPreconditionerBS(void *cvode_mem, int which, CVLsPrecSetupFnBS psetupBS, CVLsPrecSolveFnBS psolveBS)

The function CVodeSetPrecSolveFnBS() specifies the preconditioner setup and solve functions for the backward integration, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • psetupBS – user-defined preconditioner setup function.

  • psolveBS – user-defined preconditioner solve function.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

The psetupBS argument may be NULL if no setup operation is involved in the preconditioner. The previous routine CVSpilsSetPrecSolveFnBS is now deprecated.

int CVodeSetEpsLinB(void *cvode_mem, int which, sunrealtype eplifacB)

The function CVodeSetEpsLinB() specifies the factor by which the Krylov linear solver’s convergence test constant is reduced from the nonlinear iteration test constant. This routine can be used in both the cases where the backward problem does and does not depend on the forward sensitvities.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • eplifacB – value of the convergence test constant reduction factor \(\geq 0.0\).

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier, or eplifacB was negative.

Notes:

The default value is \(0.05\). Passing a value eplifacB = 0.0 also indicates using the default value. The previous routine CVSpilsSetEpsLinB is now deprecated.

int CVodeSetLSNormFactorB(void *cvode_mem, int which, sunrealtype nrmfac)

The function CVodeSetLSNormFactor() specifies the factor to use when converting from the integrator tolerance (WRMS norm) to the linear solver tolerance (L2 norm) for Newton linear system solves e.g., tol_L2 = fac * tol_WRMS. This routine can be used in both the cases wherethe backward problem does and does not depend on the forward sensitvities.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • nrmfac – the norm conversion factor. If nrmfac is: \(> 0\) then the provided value is used. \(= 0\) then the conversion factor is computed using the vector length i.e., nrmfac = N_VGetLength(y) default. \(< 0\) then the conversion factor is computed using the vector dot product nrmfac = N_VDotProd(v,v) where all the entries of v are one.

Return value:
  • CVLS_SUCCESS – The optional value has been successfully set.

  • CVLS_MEM_NULLcvode_mem was NULL.

  • CVLS_LMEM_NULL – The CVLS linear solver has not been initialized.

  • CVLS_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CVLS_ILL_INPUT – The parameter which represented an invalid identifier.

Notes:

This function must be called after the CVLS linear solver interface has been initialized through a call to CVodeSetLinearSolverB(). Prior to the introduction of N_VGetLength in SUNDIALS v5.0.0 (CVODES v5.0.0) the value of nrmfac was computed using the vector dot product i.e., the nrmfac < 0 case.

5.4.4.2.10. Optional output functions for the backward problem

The user of the adjoint module in CVODES has access to any of the optional output functions described in §5.4.1.3.12, both for the main solver and for the linear solver modules. The first argument of these CVodeGet* and CVode*Get* functions is the pointer to the CVODES memory block for the backward problem. In order to call any of these functions, the user must first call the following function to obtain this pointer.

void *CVodeGetAdjCVodeBmem(void *cvode_mem, int which)

The function CVodeGetAdjCVodeBmem() returns a pointer to the CVODES memory block for the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block created by CVodeCreate().

  • which – the identifier of the backward problem.

Return value:
  • void

Warning

The user should not modify cvode_memB in any way. Optional output calls should pass cvode_memB as the first argument; for example, to get the number of integration steps: flag = CVodeGetNumSteps(cvodes_memB, nsteps).

To get values of the forward solution during a backward integration, use the following function. The input value of t would typically be equal to that at which the backward solution has just been obtained with CVodeGetB(). In any case, it must be within the last checkpoint interval used by CVodeB().

int CVodeGetAdjY(void *cvode_mem, sunrealtype t, N_Vector y)

The function CVodeGetAdjY() returns the interpolated value of the forward solution \(y\) during a backward integration.

Arguments:
  • cvode_mem – pointer to the CVODES memory block created by CVodeCreate().

  • t – value of the independent variable at which \(y\) is desired input.

  • y – forward solution \(y(t)\).

Return value:
  • CV_SUCCESSCVodeGetAdjY() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_GETY_BADT – The value of t was outside the current checkpoint interval.

Warning

The user must allocate space for y.

int CVodeGetAdjCheckPointsInfo(void *cvode_mem, CVadjCheckPointRec *ckpnt)

The function CVodeGetAdjCheckPointsInfo() loads an array of ncheck+1 records of type CVadjCheckPointRec. The user must allocate space for the array ckpnt.

Arguments:
  • cvode_mem – pointer to the CVODES memory block created by CVodeCreate().

  • ckpnt – array of ncheck+1 checkpoint records.

Return value:
  • void

Notes:

The members of each record ckpnt[i] are:

  • ckpnt[i].my_addr (void *) – address of current checkpoint in cvode_mem->cv_adj_mem

  • ckpnt[i].next_addr (void *) – address of next checkpoint

  • ckpnt[i].t0 (sunrealtype) – start of checkpoint interval

  • ckpnt[i].t1 (sunrealtype) – end of checkpoint interval

  • ckpnt[i].nstep (long int) – step counter at ckeckpoint t0

  • ckpnt[i].order (int) – method order at checkpoint t0

  • ckpnt[i].step (sunrealtype) – step size at checkpoint t0

5.4.4.2.11. Backward integration of quadrature equations

Not only the backward problem but also the backward quadrature equations may or may not depend on the forward sensitivities. Accordingly, either CVodeQuadInitB() or CVodeQuadInitBS() should be used to allocate internal memory and to initialize backward quadratures. For any other operation (extraction, optional input/output, reinitialization, deallocation), the same function is callable regardless of whether or not the quadratures are sensitivity-dependent.

5.4.4.2.11.1. Backward quadrature initialization functions

The function CVodeQuadInitB() initializes and allocates memory for the backward integration of quadrature equations that do not depend on forward sensitivities. It has the following form:

int CVodeQuadInitB(void *cvode_mem, int which, CVQuadRhsFnB rhsQB, N_Vector yQB0)

The function CVodeQuadInitB() provides required problem specifications, allocates internal memory, and initializes backward quadrature integration.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • rhsQB – is the function which computes \(fQB\).

  • yQB0 – is the value of the quadrature variables at tB0.

Return value:
  • CV_SUCCESS – The call to CVodeQuadInitB() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_MEM_FAIL – A memory allocation request has failed.

  • CV_ILL_INPUT – The parameter which is an invalid identifier.

The function CVodeQuadInitBS() initializes and allocates memory for the backward integration of quadrature equations that depends on the forward sensitivities.

int CVodeQuadInitBS(void *cvode_mem, int which, CVQuadRhsFnBS rhsQBS, N_Vector yQBS0)

The function CVodeQuadInitBS() provides required problem specifications, allocates internal memory, and initializes backward quadrature integration.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • rhsQBS – is the function which computes \(fQBS\).

  • yQBS0 – is the value of the sensitivity-dependent quadrature variables at tB0.

Return value:
  • CV_SUCCESS – The call to CVodeQuadInitBS() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_MEM_FAIL – A memory allocation request has failed.

  • CV_ILL_INPUT – The parameter which is an invalid identifier.

The integration of quadrature equations during the backward phase can be re-initialized by calling the following function. Before calling CVodeQuadReInitB() for a new backward problem, call any desired solution extraction functions CVodeGet** associated with the previous backward problem.

int CVodeQuadReInitB(void *cvode_mem, int which, N_Vector yQB0)

The function CVodeQuadReInitB() re-initializes the backward quadrature integration.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • yQB0 – is the value of the quadrature variables at tB0.

Return value:
  • CV_SUCCESS – The call to CVodeQuadReInitB() was successful.

  • CV_MEM_NULLcvode_mem was NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_MEM_FAIL – A memory allocation request has failed.

  • CV_NO_QUAD – Quadrature integration was not activated through a previous call to CVodeQuadInitB().

  • CV_ILL_INPUT – The parameter which is an invalid identifier.

Notes:

The function CVodeQuadReInitB() can be called after a call to either CVodeQuadInitB() or CVodeQuadInitBS().

5.4.4.2.11.2. Backward quadrature extraction function

To extract the values of the quadrature variables at the last return time of CVodeB(), CVODES provides a wrapper for the function CVodeGetQuad().

int CVodeGetQuadB(void *cvode_mem, int which, sunrealtype *tret, N_Vector yQB)

The function CVodeGetQuadB() returns the quadrature solution vector after a successful return from CVodeB().

Arguments:
  • cvode_mem – pointer to the CVODES memory.

  • tret – the time reached by the solver output.

  • yQB – the computed quadrature vector.

Return value:
  • CV_SUCCESSCVodeGetQuadB() was successful.

  • CV_MEM_NULLcvode_mem is NULL.

  • CV_NO_ADJ – The function CVodeAdjInit() has not been previously called.

  • CV_NO_QUAD – Quadrature integration was not initialized.

  • CV_BAD_DKYyQB was NULL.

  • CV_ILL_INPUT – The parameter which is an invalid identifier.

Warning

The user must allocate space for yQB. To obtain the quadratures associated with a given backward problem at some other time within the last integration step, first obtain a pointer to the proper CVODES memory structure by calling CVodeGetAdjCVodeBmem() and then use it to call CVodeGetQuadDky().

5.4.4.2.11.3. Optional input/output functions for backward quadrature integration

Optional values controlling the backward integration of quadrature equations can be changed from their default values through calls to one of the following functions which are wrappers for the corresponding optional input functions defined in §5.4.2.4. The user must specify the identifier which of the backward problem for which the optional values are specified.

flag = CVodeSetQuadErrConB(cvode_mem, which, errconQ);
flag = CVodeQuadSStolerancesB(cvode_mem, which, reltolQ, abstolQ);
flag = CVodeQuadSVtolerancesB(cvode_mem, which, reltolQ, abstolQ);

Their return value flag (of type int) can have any of the return values of its counterparts, but it can also be CV_NO_ADJ if the function CVodeAdjInit() has not been previously called or CV_ILL_INPUT if the parameter which was an invalid identifier.

Access to optional outputs related to backward quadrature integration can be obtained by calling the corresponding CVodeGetQuad* functions (see §5.4.2.5). A pointer cvode_memB to the CVODES memory block for the backward problem, required as the first argument of these functions, can be obtained through a call to the functions CVodeGetAdjCVodeBmem().

5.4.4.3. User-supplied functions for adjoint sensitivity analysis

In addition to the required ODE right-hand side function and any optional functions for the forward problem, when using the adjoint sensitivity module in CVODES, the user must supply one function defining the backward problem ODE and, optionally, functions to supply Jacobian-related information and one or two functions that define the preconditioner (if an iterative SUNLinearSolver module is selected) for the backward problem. Type definitions for all these user-supplied functions are given below.

5.4.4.3.1. ODE right-hand side for the backward problem

If the backward problem does not depend on the forward sensitivities, the user must provide a rhsB function of type CVRhsFnB defined as follows:

typedef int (*CVRhsFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector yBdot, void *user_dataB)

This function evaluates the right-hand side \(f_B(t,y,y_B)\) of the backward problem ODE system. This could be either (5.21) or (5.24).

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • yBdot – is the output vector containing the right-hand side \(f_B\) of the backward ODE problem.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

A CVRhsFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_RHSFUNC_FAIL).

Notes:

Allocation of memory for yBdot is handled within CVODES. The y, yB, and yBdot arguments are all of type N_Vector, but yB and yBdot typically have different internal representations from y. It is the user’s responsibility to access the vector data consistently (including the use of the correct accessor macros from each N_Vector implementation). For the sake of computational efficiency, the vector functions in the two N_Vector implementations provided with CVODES do not perform any consistency checks with respect to their N_Vector arguments (see §9). The user_dataB pointer is passed to the user’s rhsB function every time it is called and can be the same as the user_data pointer used for the forward problem.

Warning

Before calling the user’s rhsB function, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the right-hand side function which will halt the integration and CVodeB() will return CV_RHSFUNC_FAIL.

5.4.4.3.2. ODE right-hand side for the backward problem depending on the forward sensitivities

If the backward problem does depend on the forward sensitivities, the user must provide a rhsBS function of type CVRhsFnBS defined as follows:

typedef int (*CVRhsFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector yBdot, void *user_dataB)

This function evaluates the right-hand side \(f_B(t, y, y_B, s)\) of the backward problem ODE system. This could be either (5.21) or (5.24).

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – a pointer to an array of Ns vectors containing the sensitvities of the forward solution.

  • yB – is the current value of the backward dependent variable vector.

  • yBdot – is the output vector containing the right-hand side.

  • user_dataB – is a pointer to user data, same as passed to CVodeSetUserDataB().

Return value:

A CVRhsFnBS should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_RHSFUNC_FAIL).

Notes:

Allocation of memory for qBdot is handled within CVODES. The y, yB, and yBdot arguments are all of type N_Vector, but yB and yBdot typically have different internal representations from y. Likewise for each yS[i]. It is the user’s responsibility to access the vector data consistently (including the use of the correct accessor macros from each N_Vector implementation). For the sake of computational efficiency, the vector functions in the two N_Vector implementations provided with CVODES do not perform any consistency checks with respect to their N_Vector arguments (see §9). The user_dataB pointer is passed to the user’s rhsBS function every time it is called and can be the same as the user_data pointer used for the forward problem.

Warning

Before calling the user’s rhsBS function, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the right-hand side function which will halt the integration and CVodeB() will return CV_RHSFUNC_FAIL.

5.4.4.3.3. Quadrature right-hand side for the backward problem

The user must provide an fQB function of type CVQuadRhsFnB defined by

typedef int (*CVQuadRhsFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector qBdot, void *user_dataB)

This function computes the quadrature equation right-hand side for the backward problem.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • qBdot – is the output vector containing the right-hand side fQB of the backward quadrature equations.

  • user_dataB – is a pointer to user data, same as passed to CVodeSetUserDataB().

Return value:

A CVQuadRhsFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_QRHSFUNC_FAIL).

Notes:

Allocation of memory for rhsvalBQ is handled within CVODES. The y, yB, and qBdot arguments are all of type N_Vector, but they typically do not all have the same representation. It is the user’s responsibility to access the vector data consistently (including the use of the correct accessor macros from each N_Vector implementation). For the sake of computational efficiency, the vector functions in the two N_Vector implementations provided with CVODES do not perform any consistency checks with repsect to their N_Vector arguments (see §9). The user_dataB pointer is passed to the user’s fQB function every time it is called and can be the same as the user_data pointer used for the forward problem.

Warning

Before calling the user’s fQB function, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the quadrature right-hand side function which will halt the integration and CVodeB() will return CV_QRHSFUNC_FAIL.

5.4.4.3.4. Sensitivity-dependent quadrature right-hand side for the backward problem

The user must provide an fQBS function of type CVQuadRhsFnBS defined by

typedef int (*CVQuadRhsFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector qBdot, void *user_dataB)

This function computes the quadrature equation right-hand side for the backward problem.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – a pointer to an array of Ns vectors continaing the sensitvities of the forward solution.

  • yB – is the current value of the backward dependent variable vector.

  • qBdot – is the output vector containing the right-hand side fQBS of the backward quadrature equations.

  • user_dataB – is a pointer to user data, same as passed to CVodeSetUserDataB().

Return value:

A CVQuadRhsFnBS should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_QRHSFUNC_FAIL).

Notes:

Allocation of memory for qBdot is handled within CVODES. The y, yS, and qBdot arguments are all of type N_Vector, but they typically do not all have the same internal representation. Likewise for each yS[i]. It is the user’s responsibility to access the vector data consistently (including the use of the correct accessor macros from each N_Vector implementation). For the sake of computational efficiency, the vector functions in the two N_Vector implementations provided with CVODES do not perform any consistency checks with repsect to their N_Vector arguments (see §9). The user_dataB pointer is passed to the user’s fQBS function every time it is called and can be the same as the user_data pointer used for the forward problem.

Warning

Before calling the user’s fQBS function, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the quadrature right-hand side function which will halt the integration and CVodeB() will return CV_QRHSFUNC_FAIL.

5.4.4.3.5. Jacobian construction for the backward problem (matrix-based linear solvers)

If a matrix-based linear solver module is used for the backward problem (i.e., a non-NULL SUNMatrix object was supplied to CVodeSetLinearSolverB()), the user may provide a function of type CVLsJacFnB or CVLsJacFnBS, defined as follows:

typedef int (*CVLsJacFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, SUNMatrix JacB, void *user_dataB, N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B)

This function computes the Jacobian of the backward problem (or an approximation to it).

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • JacB – is the output approximate Jacobian matrix.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmp1B, tmp2B, tmp3B – are pointers to memory allocated for variables of type N_Vector which can be used by the CVLsJacFnB function as temporary storage or work space.

Return value:

A CVLsJacFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct, while CVLS sets last_flag to CVLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

Notes:

A user-supplied Jacobian function must load the matrix JacB with an approximation to the Jacobian matrix at the point (t, y, yB), where y is the solution of the original IVP at time tt, and yB is the solution of the backward problem at the same time. 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 implementation-specific SUNMatrix interface functions (see §10 for details). With direct linear solvers (i.e., linear solvers with type SUNLINEARSOLVER_DIRECT), the Jacobian matrix \(J(t,y)\) is zeroed out prior to calling the user-supplied Jacobian function so only nonzero elements need to be loaded into JacB.

Warning

Before calling the user’s CVLsJacFnB, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the Jacobian function which will halt the integration (CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR). The previous function type CVDlsJacFnB is identical to CVLsJacFnB, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

typedef int (*CVLsJacFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, SUNMatrix JacB, void *user_dataB, N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B)

This function computes the Jacobian of the backward problem (or an approximation to it), in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – a pointer to an array of Ns vectors containing the sensitvities of the forward solution.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • JacB – is the output approximate Jacobian matrix.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmp1B, tmp2B, tmp3B – are pointers to memory allocated for variables of type N_Vector which can be used by the CVLsLinSysFnBS function as temporary storage or work space.

Return value:

A CVLsJacFnBS should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct, while CVLS sets last_flag to CVLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

Notes:

A user-supplied Jacobian function must load the matrix JacB with an approximation to the Jacobian matrix at the point (t, y, yS, yB), where y is the solution of the original IVP at time tt, yS is the vector of forward sensitivities at time tt, and yB is the solution of the backward problem at the same time. 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 implementation-specific SUNMatrix interface functions (see §10). With direct linear solvers (i.e., linear solvers with type SUNLINEARSOLVER_DIRECT, the Jacobian matrix \(J(t,y)\) is zeroed out prior to calling the user-supplied Jacobian function so only nonzero elements need to be loaded into JacB.

Warning

Before calling the user’s CVLsJacFnBS, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the Jacobian function which will halt the integration (CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR). The previous function type CVDlsJacFnBS is identical to CVLsJacFnBS, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

5.4.4.3.6. Linear system construction for the backward problem (matrix-based linear solvers)

With matrix-based linear solver modules, as an alternative to optionally supplying a function for evaluating the Jacobian of the ODE right-hand side function, the user may optionally supply a function of type CVLsLinSysFnB or CVLsLinSysFnBS for evaluating the linear system, \(M_B = I - \gamma_B J_B\) (or an approximation of it) for the backward problem.

typedef int (*CVLsLinSysFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, SUNMatrix AB, sunbooleantype jokB, sunbooleantype *jcurB, sunrealtype gammaB, void *user_dataB, N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B);

This function computes the linear system of the backward problem (or an approximation to it).

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • AB – is the output approximate linear system matrix.

  • jokB – is an input flag indicating whether Jacobian-related data needs to be recomputed (jokB = SUNFALSE) or informtion saved from a previous information can be safely used (jokB = SUNTRUE).

  • jcurB – is an output flag which must be set to SUNTRUE if Jacobian-related data was recomputed or SUNFALSE otherwise.

  • gammaB – is the scalar appearing in the matrix \(M_B = I - \gamma_B J_B\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmp1B, tmp2B, tmp3B – are pointers to memory allocated for variables of type N_Vector which can be used by the CVLsLinSysFnB function as temporary storage or work space.

Return value:

A CVLsLinSysFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct, while CVLS sets last_flag to CVLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

Notes:

A user-supplied linear system function must load the matrix AB with an approximation to the linear system matrix at the point (t, y, yB), where y is the solution of the original IVP at time tt, and yB is the solution of the backward problem at the same time.

Warning

Before calling the user’s CVLsLinSysFnB, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the linear system function which will halt the integration (CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

typedef int (*CVLsLinSysFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, SUNMatrix AB, sunbooleantype jokB, sunbooleantype *jcurB, sunrealtype gammaB, void *user_dataB, N_Vector tmp1B, N_Vector tmp2B, N_Vector tmp3B);

This function computes the linear system of the backward problem (or an approximation to it), in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – a pointer to an array of Ns vectors containing the sensitivities of the forward solution.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • AB – is the output approximate linear system matrix.

  • jokB – is an input flag indicating whether Jacobian-related data needs to be recomputed (jokB = SUNFALSE) or informtion saved from a previous information can be safely used (jokB = SUNTRUE).

  • jcurB – is an output flag which must be set to SUNTRUE if Jacobian-related data was recomputed or SUNFALSE otherwise.

  • gammaB – is the scalar appearing in the matrix

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmp1B, tmp2B, tmp3B – are pointers to memory allocated for variables of type N_Vector which can be used by the CVLsLinSysFnBS function as temporary storage or work space.

Return value:

A CVLsLinSysFnBS should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct, while CVLS sets last_flag to CVLS_JACFUNC_RECVR), or a negative value if it failed unrecoverably (in which case the integration is halted, CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

Notes:

A user-supplied linear system function must load the matrix AB with an approximation to the linear system matrix at the point (t, y, yS, yB), where y is the solution of the original IVP at time tt, yS is the vector of forward sensitivities at time t, and yB is the solution of the backward problem at the same time.

Warning

Before calling the user’s CVLsLinSysFnBS, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the linear system function which will halt the integration (CVodeB() returns CV_LSETUP_FAIL and CVLS sets last_flag to CVLS_JACFUNC_UNRECVR).

5.4.4.3.7. Jacobian-vector product for the backward problem (matrix-free linear solvers)

If a matrix-free linear solver is to be used for the backward problem (i.e., a NULL-valued SUNMatrix was supplied to CVodeSetLinearSolverB() in the steps described in §5.4.4.1, the user may provide a function of type CVLsJacTimesVecFnB or CVLsJacTimesVecFnBS in the following form, to compute matrix-vector products \(Jv\). If such a function is not supplied, the default is a difference quotient approximation to these products.

typedef int (*CVLsJacTimesVecFnB)(N_Vector vB, N_Vector JvB, sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, void *jac_dataB, N_Vector tmpB);

This function computes the action of the Jacobian JB for the backward problem on a given vector vB.

Arguments:
  • vB – is the vector by which the Jacobian must be multiplied to the right.

  • JvB – is the computed output vector JB*vB.

  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmpB – is a pointer to memory allocated for a variable of type N_Vector which can be used by CVLsJacTimesVecFnB as temporary storage or work space.

Return value:

The return value of a function of type CVLsJacTimesVecFnB should be if successful or nonzero if an error was encountered, in which case the integration is halted.

Notes:

A user-supplied Jacobian-vector product function must load the vector JvB with the product of the Jacobian of the backward problem at the point (t, y, yB) and the vector vB. Here, y is the solution of the original IVP at time t and yB is the solution of the backward problem at the same time. The rest of the arguments are equivalent to those passed to a function of type CVLsJacTimesVecFn. If the backward problem is the adjoint of \({\dot y} = f(t, y)\), then this function is to compute \(-({\partial f}/{\partial y_i})^T v_B\). The previous function type CVSpilsJacTimesVecFnB is deprecated.

typedef int (*CVLsJacTimesVecFnBS)(N_Vector vB, N_Vector JvB, sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, void *user_dataB, N_Vector tmpB);

This function computes the action of the Jacobian JB for the backward problem on a given vector vB, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • vB – is the vector by which the Jacobian must be multiplied to the right.

  • JvB – is the computed output vector JB*vB.

  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – is a pointer to an array containing the forward sensitivity vectors.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

  • tmpB – is a pointer to memory allocated for a variable of type N_Vector which can be used by CVLsJacTimesVecFnB as temporary storage or work space.

Return value:

The return value of a function of type CVLsJacTimesVecFnBS should be if successful or nonzero if an error was encountered, in which case the integration is halted.

Notes:

A user-supplied Jacobian-vector product function must load the vector JvB with the product of the Jacobian of the backward problem at the point (t, y, yB) and the vector vB. Here, y is the solution of the original IVP at time t and yB is the solution of the backward problem at the same time. The rest of the arguments are equivalent to those passed to a function of type CVLsJacTimesVecFn. The previous function type CVSpilsJacTimesVecFnBS is deprecated.

5.4.4.3.8. Jacobian-vector product setup for the backward problem (matrix-free linear solvers)

If the user’s Jacobian-times-vector routine requires that any Jacobian-related data be preprocessed or evaluated, then this needs to be done in a user-supplied function of type CVLsJacTimesSetupFnB or CVLsJacTimesSetupFnBS, defined as follows:

typedef int (*CVLsJacTimesSetupFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, void *user_dataB)

This function preprocesses and/or evaluates Jacobian data needed by the Jacobian-times-vector routine for the backward problem.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the dependent variable vector, \(y(t)\).

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the right-hand-side for the backward problem.

  • user_dataB – is a pointer to user data CVodeSetUserDataB().

Return value:

The value returned by the Jacobian-vector setup function should be if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes:

Each call to the Jacobian-vector setup function is preceded by a call to the backward problem residual user function with the same (t,y, yB) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the right-hand-side function. If the user’s CVLsJacTimesVecFnB function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to cvode_mem to user_dataB and then use the CVGet* functions described in §5.4.1.3.12. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h. The previous function type CVSpilsJacTimesSetupFnB is identical to CVLsJacTimesSetupFnB, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

typedef int (*CVLsJacTimesSetupFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, void *user_dataB)

This function preprocesses and/or evaluates Jacobian data needed by the Jacobian-times-vector routine for the backward problem, in the case that the backward problem depends on the forward sensitivities.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the dependent variable vector, \(y(t)\).

  • yS – a pointer to an array of Ns vectors containing the sensitvities of the forward solution.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the right-hand-side function for the backward problem.

  • user_dataB – is a pointer to the same user data provided to CVodeSetUserDataB().

Return value:

The value returned by the Jacobian-vector setup function should be if successful, positive for a recoverable error (in which case the step will be retried), or negative for an unrecoverable error (in which case the integration is halted).

Notes:

Each call to the Jacobian-vector setup function is preceded by a call to the backward problem residual user function with the same (t,y, yS, yB) arguments. Thus, the setup function can use any auxiliary data that is computed and saved during the evaluation of the right-hand-side function. If the user’s CVLsJacTimesVecFnBS function uses difference quotient approximations, it may need to access quantities not in the call list. These include the current stepsize, the error weights, etc. To obtain these, the user will need to add a pointer to cvode_mem to user_dataB and then use the CVGet* functions described in §5.4.1.3.12. The unit roundoff can be accessed as SUN_UNIT_ROUNDOFF defined in sundials_types.h. The previous function type CVSpilsJacTimesSetupFnBS is identical to CVLsJacTimesSetupFnBS, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

5.4.4.3.9. Preconditioner solve for the backward problem (iterative linear solvers)

If a user-supplied 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\) may be either a left or a right preconditioner matrix. Here \(P\) should approximate (at least crudely) the matrix \(M_B = I - \gamma_B J_B\), where \(J_B = \partial f_B/ \partial y_B\). If preconditioning is done on both sides, the product of the two preconditioner matrices should approximate \(M_B\). This function must be of one of the following two types:

typedef int (*CVLsPrecSolveFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, N_Vector rvecB, N_Vector zvecB, sunrealtype gammaB, sunrealtype deltaB, void *user_dataB)

This function solves the preconditioning system \(Pz = r\) for the backward problem.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • rvecB – is the right-hand side vector r of the linear system to be solved.

  • zvecB – is the computed output vector.

  • gammaB – is the scalar appearing in the matrix, \(M_B = I - \gamma_B J_B\).

  • deltaB – is an input tolerance to be used if an iterative method is employed in the solution.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

The return value of a preconditioner solve function for the backward problem should be 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 previous function type CVSpilsPrecSolveFnB is deprecated.

typedef int (*CVLsPrecSolveFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, N_Vector rvecB, N_Vector zvecB, sunrealtype gammaB, sunrealtype deltaB, void *user_dataB)

This function solves the preconditioning system \(Pz = r\) for the backward problem, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – is a pointer to an array containing the forward sensitivity vectors.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • rvecB – is the right-hand side vector r of the linear system to be solved.

  • zvecB – is the computed output vector.

  • gammaB – is the scalar appearing in the matrix, \(M_B = I - \gamma_B J_B\).

  • deltaB – is an input tolerance to be used if an iterative method is employed in the solution.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

The return value of a preconditioner solve function for the backward problem should be 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 previous function type CVSpilsPrecSolveFnBS is identical to CVLsPrecSolveFnBS, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

5.4.4.3.10. Preconditioner setup for the backward problem (iterative linear solvers)

If the user’s preconditioner requires that any Jacobian-related data be preprocessed or evaluated, then this needs to be done in a user-supplied function of one of the following two types:

typedef int (*CVLsPrecSetupFnB)(sunrealtype t, N_Vector y, N_Vector yB, N_Vector fyB, sunbooleantype jokB, sunbooleantype *jcurPtrB, sunrealtype gammaB, void *user_dataB)

This function preprocesses and/or evaluates Jacobian-related data needed by the preconditioner for the backward problem.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • jokB – is an input flag indicating whether Jacobian-related data needs to be recomputed (jokB = SUNFALSE) or information saved from a previous invokation can be safely used (jokB = SUNTRUE).

  • jcurPtr – is an output flag which must be set to SUNTRUE if Jacobian-related data was recomputed or SUNFALSE otherwise.

  • gammaB – is the scalar appearing in the matrix \(M_B = I - \gamma_B J_B\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

The return value of a preconditioner setup function for the backward problem should be 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 previous function type CVSpilsPrecSetupFnB is identical to CVLsPrecSetupFnB, and may still be used for backward-compatibility. However, this will be deprecated in future releases, so we recommend that users transition to the new function type name soon.

typedef int (*CVLsPrecSetupFnBS)(sunrealtype t, N_Vector y, N_Vector *yS, N_Vector yB, N_Vector fyB, sunbooleantype jokB, sunbooleantype *jcurPtrB, sunrealtype gammaB, void *user_dataB)

This function preprocesses and/or evaluates Jacobian-related data needed by the preconditioner for the backward problem, in the case where the backward problem depends on the forward sensitivities.

Arguments:
  • t – is the current value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yS – is a pointer to an array containing the forward sensitivity vectors.

  • yB – is the current value of the backward dependent variable vector.

  • fyB – is the current value of the backward right-hand side function \(f_B\).

  • jokB – is an input flag indicating whether Jacobian-related data needs to be recomputed (jokB = SUNFALSE) or information saved from a previous invokation can be safely used (jokB = SUNTRUE).

  • jcurPtr – is an output flag which must be set to SUNTRUE if Jacobian-related data was recomputed or SUNFALSE otherwise.

  • gammaB – is the scalar appearing in the matrix \(M_B = I - \gamma_B J_B\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

The return value of a preconditioner setup function for the backward problem should be 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 previous function type CVSpilsPrecSetupFnBS is deprecated.

5.4.4.4. Using CVODES preconditioner modules for the backward problem

As on the forward integration phase, the efficiency of Krylov iterative methods for the solution of linear systems can be greatly enhanced through preconditioning. Both preconditioner modules provided with SUNDIALS, the serial banded preconditioner CVBANDPRE and the parallel band-block-diagonal preconditioner module CVBBDPRE, provide interface functions through which they can be used on the backward integration phase.

5.4.4.4.1. Using the banded preconditioner CVBANDPRE

The adjoint module in CVODES offers an interface to the banded preconditioner module CVBANDPRE described in section §5.4.2.7.1. This preconditioner, usable only in a serial setting, provides a band matrix preconditioner based on difference quotients of the backward problem right-hand side function fB. It generates a banded approximation to the Jacobian with \(m_{lB}\) sub-diagonals and \(m_{uB}\) super-diagonals to be used with one of the Krylov linear solvers.

In order to use the CVBANDPRE module in the solution of the backward problem, the user need not define any additional functions. Instead, after an iterative SUNLinearSolver object has been attached to CVODES via a call to CVodeSetLinearSolverB(), the following call to the CVBANDPRE module initialization function must be made.

int CVBandPrecInitB(void *cvode_mem, int which, sunindextype nB, sunindextype muB, sunindextype mlB)

The function CVBandPrecInitB() initializes and allocates memory for the CVBANDPRE preconditioner for the backward problem. It creates, allocates, and stores (internally in the CVODES solver block) a pointer to the newly created CVBANDPRE memory block.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • nB – backward problem dimension.

  • muB – upper half-bandwidth of the backward problem Jacobian approximation.

  • mlB – lower half-bandwidth of the backward problem Jacobian approximation.

Return value:
  • CVLS_SUCCESS – The call to CVodeBandPrecInitB() was successful.

  • CVLS_MEM_FAIL – A memory allocation request has failed.

  • CVLS_MEM_NULL – The cvode_mem argument was NULL.

  • CVLS_LMEM_NULL – No linear solver has been attached.

  • CVLS_ILL_INPUT – An invalid parameter has been passed.

For more details on CVBANDPRE see §5.4.2.7.1.

5.4.4.4.2. Using the band-block-diagonal preconditioner CVBBDPRE

The adjoint module in CVODES offers an interface to the band-block-diagonal preconditioner module CVBBDPRE described in section §5.4.2.7.2. This generates a preconditioner that is a block-diagonal matrix with each block being a band matrix and can be used with one of the Krylov linear solvers and with the MPI-parallel vector module NVECTOR_PARALLEL.

In order to use the CVBBDPRE module in the solution of the backward problem, the user must define one or two additional functions, described at the end of this section.

5.4.4.4.2.1. Initialization of CVBBDPRE

The CVBBDPRE module is initialized by calling the following function, after an iterative SUNLinearSolver object has been attached to CVODES via a call to CVodeSetLinearSolverB().

int CVBBDPrecInitB(void *cvode_mem, int which, sunindextype NlocalB, sunindextype mudqB, sunindextype mldqB, sunindextype mukeepB, sunindextype mlkeepB, sunrealtype dqrelyB, CVBBDLocalFnB glocB, CVBBDCommFnB gcommB)

The function CVBBDPrecInitB() initializes and allocates memory for the CVBBDPRE preconditioner for the backward problem. It creates, allocates, and stores (internally in the CVODES solver block) a pointer to the newly created CVBBDPRE memory block.

Arguments:
  • cvode_mem – pointer to the CVODES memory block.

  • which – the identifier of the backward problem.

  • NlocalB – local vector dimension for the backward problem.

  • mudqB – upper half-bandwidth to be used in the difference-quotient Jacobian approximation.

  • mldqB – lower half-bandwidth to be used in the difference-quotient Jacobian approximation.

  • mukeepB – upper half-bandwidth of the retained banded approximate Jacobian block.

  • mlkeepB – lower half-bandwidth of the retained banded approximate Jacobian block.

  • dqrelyB – the relative increment in components of yB used in the difference quotient approximations. The default is \(\text{dqrelyB} = \sqrt{\text{unit roundoff}}\) , which can be specified by passing dqrely = 0.0.

  • glocB – the function which computes the function \(g_Bt,y,y_B\) approximating the right-hand side of the backward problem.

  • gcommB – the optional function which performs all interprocess communication required for the computation of \(g_B\).

Return value:
  • CVLS_SUCCESS – The call to CVodeBBDPrecInitB() was successful.

  • CVLS_MEM_FAIL – A memory allocation request has failed.

  • CVLS_MEM_NULL – The cvode_mem argument was NULL.

  • CVLS_LMEM_NULL – No linear solver has been attached.

  • CVLS_ILL_INPUT – An invalid parameter has been passed.

int CVBBDPrecReInitB(void *cvode_mem, int which, sunindextype mudqB, sunindextype mldqB, sunrealtype dqrelyB)

The function CVBBDPrecReInitB() reinitializes the CVBBDPRE preconditioner for the backward problem.

Arguments:
  • cvode_mem – pointer to the CVODES memory block returned by CVodeCreate().

  • which – the identifier of the backward problem.

  • mudqB – upper half-bandwidth to be used in the difference-quotient Jacobian approximation.

  • mldqB – lower half-bandwidth to be used in the difference-quotient Jacobian approximation.

  • dqrelyB – the relative increment in components of yB used in the difference quotient approximations.

Return value:
  • CVLS_SUCCESS – The call to CVodeBBDPrecReInitB() was successful.

  • CVLS_MEM_FAIL – A memory allocation request has failed.

  • CVLS_MEM_NULL – The cvode_mem argument was NULL.

  • CVLS_PMEM_NULL – The CVodeBBDPrecInitB() has not been previously called.

  • CVLS_LMEM_NULL – No linear solver has been attached.

  • CVLS_ILL_INPUT – An invalid parameter has been passed.

For more details on CVBBDPRE see §5.4.2.7.2.

5.4.4.4.2.2. User-supplied functions for CVBBDPRE

To use the CVBBDPRE module, the user must supply one or two functions which the module calls to construct the preconditioner: a required function glocB (of type CVBBDLocalFnB) which approximates the right-hand side of the backward problem and which is computed locally, and an optional function gcommB (of type CVBBDCommFnB) which performs all interprocess communication necessary to evaluate this approximate right-hand side. The prototypes for these two functions are described below.

typedef int (*CVBBDLocalFnB)(sunindextype NlocalB, sunrealtype t, N_Vector y, N_Vector yB, N_Vector gB, void *user_dataB)

This glocB function loads the vector gB, an approximation to the right-hand side \(f_B\) of the backward problem, as a function of t, y, and yB.

Arguments:
  • NlocalB – is the local vector length for the backward problem.

  • t – is the value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • gB – is the output vector, \(g_B(t, y, y_B)\).

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

An CVBBDLocalFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_LSETUP_FAIL).

Notes:

This routine must assume that all interprocess communication of data needed to calculate gB has already been done, and this data is accessible within user_dataB.

Warning

Before calling the user’s CVBBDLocalFnB, CVODES needs to evaluate (through interpolation) the values of the states from the forward integration. If an error occurs in the interpolation, CVODES triggers an unrecoverable failure in the preconditioner setup function which will halt the integration (CVodeB() returns CV_LSETUP_FAIL).

typedef int (*CVBBDCommFnB)(sunindextype NlocalB, sunrealtype t, N_Vector y, N_Vector yB, void *user_dataB)

This gcommB function must perform all interprocess communications necessary for the execution of the glocB function above, using the input vectors y and yB.

Arguments:
  • NlocalB – is the local vector length.

  • t – is the value of the independent variable.

  • y – is the current value of the forward solution vector.

  • yB – is the current value of the backward dependent variable vector.

  • user_dataB – is a pointer to the same user data passed to CVodeSetUserDataB().

Return value:

An CVBBDCommFnB should return 0 if successful, a positive value if a recoverable error occurred (in which case CVODES will attempt to correct), or a negative value if it failed unrecoverably (in which case the integration is halted and CVodeB() returns CV_LSETUP_FAIL).

Notes:

The gcommB function is expected to save communicated data in space defined within the structure user_dataB. Each call to the gcommB function is preceded by a call to the function that evaluates the right-hand side of the backward problem with the same t, y, and yB, arguments. If there is no additional communication needed, then pass gcommB = NULL to CVBBDPrecInitB().