cvode(),
cvodes(), ida() and cvsolve()
accept an optional jacobian argument, an R
function giving the Jacobian of the system analytically. When it is not
supplied, which remains the default, SUNDIALS approximates
the Jacobian by finite differences as before, so existing code is
unaffected. For cvode(), cvodes() and
cvsolve() the function has the same signature as the system
itself, function(t, y, p), and returns an n-by-n matrix
whose [i, j] entry is d(ydot[i])/d(y[j]). For
ida() it is function(t, y, ydot, cj, p) and
returns dF/dy + cj * dF/dydot, cj being a
scalar supplied by the solver. Supplying the Jacobian is usually faster
and more accurate on stiff systemsinst/include/sundialr_capi.h and registered as
R_RegisterCCallable entry points under the package name
sundialr. A consumer reaches it with
R_GetCCallable("sundialr", ...) after
Imports: sundialr; the header exposes only
double/int and an opaque handle, so no
SUNDIALS headers or LinkingTo are required. The handle owns
a persistent state vector and is reused across many segments via
sundialr_cvode_reinit(), the pattern a host integrator loop
needs; tolerances may be scalar or per-equation (the latter realised
with a custom error-weight function, so vector relative tolerances are
supported), an analytic Jacobian is optional, and every entry point
returns a CVODE status code and records a retrievable
message rather than throwing across the foreign call boundary. It lets
CVODE be embedded as an integrator inside another package’s
own solve loop; the existing R interface is unchangedcvodes() accepts an
optional sensitivity argument, an R function
giving the sensitivity right-hand side analytically. When it is not
supplied, which remains the default, SUNDIALS approximates
the sensitivity equations by finite differences of the state right-hand
side as before, so existing code is unaffected. The function has
signature function(t, y, ydot, iS, yS, p) and returns
d(yS_iS)/dt = J %*% yS_iS + df/dp_iS as a numeric vector of
length(y), where iS is the 1-based index of
the parameter. It is independent of the jacobian argument —
the caller embeds the Jacobian inside their own sensitivity
function — so supplying one does not require supplying the other.
Providing the sensitivity right-hand side is usually faster and more
accurate than the finite-difference default, which evaluates the state
right-hand side an extra length(p) times per internal
stepsrc/scripts/cran_patches.sh) which now verifies that no
flagged call survives patchingtools/strip_sundials_tarball.sh to reproducibly
strip the bundled SUNDIALS tarball of examples, docs, tests and
benchmarksEvents record at
time 0 in cvsolve() now adds its value to the initial
condition instead of replacing it, so the IC supplied by
the caller is always the starting state and t = 0 behaves like every
other event time. Several events at t = 0 on one state now accumulate
rather than the last one winning. Results change for any call that doses
at time 0: with IC = c(1) and a value of 100 at t = 0 the
state now starts at 101 rather than 100, and
inst/examples/cvsolve_1D.r starts at 101 accordingly. To
reproduce the old numbers, subtract the initial condition from the value
of any event at time 0cvsolve() treating the initial time as 0 rather
than as the first element of time_vector. An event at the
initial time is now applied to the initial conditions as intended;
previously, if the time vector began anywhere other than 0, such an
event was silently dropped and its value lost with no error or
warningcvsolve() now rejects an event time outside the range
of time_vector instead of accepting it. An event before the
first output time cannot be applied, since the solve starts there; it
previously produced a row of zeros in the returned matrix and discarded
the event. Event times that are NA are rejected for the
same reasoncvsolve() no longer returns rows of zeros. A record
that does not advance the solution — one at the initial time, or one
sharing a time with the record before it — was skipped without storing
anything, leaving that row filled with zeros instead of the stateida() and cvsolve() where a
scalar abstolerance (the default) allocated a length-1
tolerance vector instead of one entry per state. On systems with more
than one state this read and wrote past the end of that vector, and the
solvers failed immediately with a corrector convergence error at the
initial time. Passing abstolerance as a vector of the same
length as IC was unaffected, as were cvode()
and cvodes()cvode(), cvodes(), ida() and
cvsolve() no longer leak memory when a call ends in an
error. The SUNDIALS objects behind each solve are now released by a
scope guard instead of by free calls at the end of the function, so they
are also freed when an error unwinds the stack. Previously a failed call
leaked everything allocated up to the point of the error, which
accumulated when solvers were called repeatedly (for example from an
optimiser or a fitting loop). This includes errors raised by SUNDIALS
itself, such as a failure to converge: these previously aborted the call
in a way that skipped all cleanup, and were the most common source of
the leak, since a solver that fails to converge is exactly the case a
user is likely to retry in a loop. Separately, cvsolve()
never freed its constraints vector on any path, including successful
onescvsolve() now rejects an invalid state index in the
first column of Events. The index is converted to a
position in the state vector and used to subscript it directly, so a
value outside 1:length(IC) wrote outside that vector: an
index of 0, the natural mistake when counting from zero,
corrupted the heap and could abort the R session.
Fractional and NA indices are rejected for the same reason,
a fractional index having previously been truncated to a different state
without warningR function is
now returned to the solver as a failure code rather than being thrown
out through SUNDIALS’ own C code, which lets SUNDIALS unwind its
internal state instead of being abandoned part-way through a step. The
R error itself is unchanged: the original condition still
reaches the caller, so tryCatch() on a specific condition
class continues to workR function returning the wrong number
of values is now reported instead of being read past the end. Previously
a right-hand side or residual function returning fewer values than there
are states, or a Jacobian function returning a matrix of the wrong
shape, was copied out element by element regardless and the solve
continued on uninitialised memoryida() and
cvsolve(), which previously had none, including a check
that a scalar and an equivalent vector abstolerance give
identical results, and for the callback error handling aboveCVodeSetLinearSolver was reported
as CVDlsSetLinearSolver (a name removed in SUNDIALS 4) in
cvode() and cvsolve(),
CVodeWFtolerances as CVodeSetEwtFn and
CVodeSensInit1 as CVodeSensInit in
cvodes(), and IDASetUserData and
IDASolve as their cvode() equivalents in
ida()/tools/cmake_call.sh file to remove warnings
with respect to stdout/stderr/fprintf etc in upstream C libraryUpdated build system to use cmake
Bug fixes: * Fixed the bug in assigning absolute tolerances
Bug fixes: * Fixed the linking problem due to multiple defined symbols
New features: * added the cvsolve function to solve ODEs
with multiple discontinuities
Bug fixes: * None
API change: * NONE is existing functions * A new function
cvsolve is added
Internal changes: * Now uses SUNDIALS version 5.2.0 (released March 2020) at the back end