YABI — Yet Another BDF Integrator

Note on origin: This is a complete rewrite of the BDF integrator developed by Dominik Skanda (based on his dissertation work at University of Freiburg). The mathematical formulation follows BDF methods with sensitivity generation as presented in the dissertation by Dominik Skanda.
It supports dense (LAPACK), sparse (UMFPACK/KLU via SuiteSparse), and matrix-free Newton–Krylov (GMRES/FGMRES/BiCGSTAB via Jacobian–vector products) linear solvers through a pluggable LinearSolver interface. The rewrite was carried out using OpenCode with the Big Pickle model and OpenRouter free models. All algorithmic components (Nordsieck predictor, simplified Newton, matrix scaling/conditioning strategy, corrector acceptance, and sensitivity replay) were reimplemented to match the mathematical structure of the original approach while exposing a clean modern interface.

YABI is a variable-order, variable-step BDF integrator for stiff ordinary differential equations. It provides forward, adjoint (reverse) and simultaneous sensitivity analysis through CppAD automatic differentiation.

• Orders 1–6, variable order and variable step size.
• Three right-hand-side forms: autonomous f(y), time dependent f(t, y), and time/parameter dependent f(t, y, p).
• Sensitivities: forward mode, adjoint/reverse mode, and simultaneous (variational) sensitivities.
• Stiff problems: implicit corrector with Newton iteration, Jacobian reuse, and LAPACK linear solves.
• Python bindings: import yabi — zero-copy NumPy I/O with an analytic Jacobian (v1, analytic path; CppAD AD and sensitivities in Python are the next milestone — see "Python bindings" below).
• C++17, BLAS/LAPACK, and an external CppAD installation.

yabi-torch — pip install yabi-torch

A batched, differentiable solve_ivp for torch, on top of the BDF + continuous adjoint:

pip install yabi-torch

import torch
from yabi_torch import solve_ivp

def f(y):  # y' = f(y) — any torch.nn.Module or callable
    return torch.stack([-0.04 * y[0] + 1e4 * y[1] * y[2],
                        0.04 * y[0] - 1e4 * y[1] * y[2] - 3e7 * y[1] ** 2,
                        3e7 * y[1] ** 2])

y0 = torch.tensor([1.0, 0.0, 0.0], dtype=torch.float64, requires_grad=True)
yT = solve_ivp(f, y0, t_span=(0.0, 1e5), rtol=1e-6, atol=1e-10)
yT.sum().backward()   # y0.grad = d(sum y(T)) / d y0 via the continuous adjoint

f is a torch.nn.Module (or any callable y -> f(y) on a float64 tensor); the Jacobian is automatic (torch.autograd) unless you pass an analytic jac. See "Python bindings" below for the full API (solver="dense"|"matrix_free", threads=N for a parallel batch, device="cuda").

Quick Verification

Prove tolerance adherence, first-order convergence, and higher-order sensitivity in one command. This regenerates the work–precision diagram (work_precision_diagram.png) and prints the slope summary plus the higher-order (forward/reverse, orders 1–4) sensitivity output locally.

cmake -S . -B build && cmake --build build -j
python3 reproduce_benchmarks.py

reproduce_benchmarks.py runs two self-checking testbench binaries and parses their output:

A broader, machine-readable performance/accuracy table (runtime, Jacobian evaluations, factorizations, steps, order, and error vs an independent SciPy Radau reference) is generated by the telemetry harness:

cmake --build build --target yabi_telemetry
python3 benchmarks/report.py build/yabi_telemetry

This writes benchmarks/TELEMETRY.md (the table) and benchmarks/work_precision.png (the error-vs-rtol plot). See benchmarks/ and the reference values recorded in benchmarks/telemetry.cpp.

• build/yabi_test_11_tolerance_and_order — the work–precision proof (tolerance adherence + first-order step convergence).
• build/yabi_test_28_second_order_sensitivity — exact forward and reverse sensitivities of orders 1–4.

Plotting requires matplotlib (the build and the two binaries themselves do not):

pip install matplotlib

Dependencies

Dependency	Notes
CppAD	Optional. For automatic differentiation / examples; can be disabled or not required in some builds. Must provide ADFun::capacity_order (2013 or newer) if used.
BLAS	Fortran dgemm_, dgemv_, dger_, daxpy_, dscal_, dnrm2_.
LAPACK	Fortran dgetrf_, dgetrs_ (for dense linear solver).
SuiteSparse (UMFPACK/KLU)	Optional, for sparse linear solvers via the pluggable LinearSolver interface.
C++17 compiler	GCC or Clang.
CMake ≥ 3.16	Build only.
pybind11 ≥ 3.0	Optional, for the yabi Python module (-DYABI_BUILD_PYTHON=ON). pip install pybind11.
python3-dev	Optional, for the Python module (Python.h headers).
matplotlib	Optional. Only for reproduce_benchmarks.py's work–precision plot (and pyplot/plot_all.py).

CppAD is distributed under EPL-2.0 with a GPL-2.0-or-later secondary license. BLAS and LAPACK are BSD-style. YABI itself is Apache-2.0.

On Debian/Ubuntu:

sudo apt install build-essential liblapack-dev libblas-dev cmake libcppad-dev
# optional, for the Python bindings (T37):
sudo apt install python3-dev && pip install pybind11

Layout

yabi-bdf/
  include/yabi/
    bdf_integrator.hpp        public class declaration
    linear_algebra.hpp        BLAS/LAPACK bindings
    detail/
      core.ipp                construction, reset, configuration, results
      step_control.ipp        step size, predictor, order control, output
      newton.ipp              Newton iteration and corrector acceptance
      sensitivity_setup.ipp   sensitivity storage and tape lifetime
      sensitivity_forward.ipp forward-mode replay
      sensitivity_reverse.ipp reverse/adjoint replay
  src/
    bdf_integrator.cpp        explicit instantiation for double
  examples/
  tests/
    unit/                     isolated component unit tests
    component/                cross-component integration tests
    system/                   mathematical benchmarks and regressions

The .ipp files contain the member definitions of the class template. They are implementation details: include <yabi/bdf_integrator.hpp>, not the .ipp files.

Build

cmake -S . -B build [-DCPPAD_ROOT=/path/to/cppad]
cmake --build build -j
ctest --test-dir build --output-on-failure

If CppAD is not installed in standard paths:

cmake -S . -B build -DCPPAD_ROOT=/path/to/cppad
cmake --build build -j$(nproc)
ctest --test-dir build

When CppAD is installed system-wide, -DCPPAD_ROOT / CPPAD_ROOT can be omitted. CMake first tries pkg-config, then looks for cppad/cppad.hpp and libcppad_lib.

Python bindings (T37)

An optional yabi Python module exposes the integrator on the analytic path: the RHS and analytic Jacobian are Python callables, and the state buffers are zero-copy NumPy arrays.

pip install pybind11         # >= 3.0 (Python 3.14)
sudo apt install python3-dev # Python.h headers
cmake -S . -B build-py -DYABI_BUILD_PYTHON=ON
cmake --build build-py --target yabi_python
ctest --test-dir build-py -R python_bindings

import numpy as np, yabi
def f(y):   return np.array([...])   # RHS, n doubles
def jac(y): return np.array([...])   # column-major J, n*n doubles
it = yabi.BDFIntegrator(f, jac, n)
y0 = np.array([...]); atol = np.array([...])
it.reset_integration(y0, max_step, rtol, atol)   # zero-copy y0
it.set_algorithm_settings(5, np.finfo(float).eps, True, 1)
it.integrate(t_end)
it.get_solution(y_out)                            # zero-copy write into y_out

bindings/python/test_yabi.py is the Brusselator smoke test; bindings/python/test_autograd.py checks the continuous-adjoint binding against central differences.

The first-order continuous adjoint (stable reverse sensitivity) is exposed as integrator.continuous_adjoint(seed_weights, adj_y0) — the backward hook for autograd wrappers. bindings/python/torch_wrapper.py wraps it in a torch.autograd.Function (forward solve + adjoint backward), so gradients flow through the ODE solve w.r.t. the initial condition. bindings/python/yabi_torch.py provides a batched, differentiable solve_ivp(f, y0, t_span, rtol, atol, jac=None, solver="dense", threads=1, device="cpu") on top of it: f is a torch.nn.Module (or any torch callable y -> f(y)), the Jacobian is automatic from torch.autograd unless an analytic jac is given, and solver="matrix_free" drives the Newton solve with a torch.autograd JVP (forward-mode, one J·v per Krylov iteration) instead of the full n×n Jacobian. y0 may be a batch; the gradient is the continuous adjoint. threads=N solves the batch's independent trajectories concurrently (the binding releases the GIL during the C++ solve and re-acquires it only for the Python callbacks), bit-identical to serial. device="cuda" offloads the dense Newton LU to cuSOLVER (CudaLinearSolver) and, when the GPU has compute capability ≥ 7.5, evaluates f/the Jacobian on-device (torch CUDA); it requires a CUDA build of the module (YABI_USE_CUDA=ON). Honest scope: the BDF state stays host-side (host↔device copy per evaluation); a full device-resident-state BDF is a larger follow-up. The cuSOLVER LU offload is verified bit-identical to CPU; the f-on-GPU path is implemented but unverified on the development GPU (a Pascal GTX 1050, CC 6.1, which torch's cu12x wheels no longer support). test_solve_ivp.py runs the torch.autograd.gradcheck, batch (forward + adjoint), and stiff Robertson checks (dense + matrix-free); test_nn_module.py checks an nn.Module f against the matrix exponential and cross-mode agreement; test_batch_parallel.py checks threads=N is bit-identical to serial; test_cuda.py checks the cuSOLVER LU offload and (gated on CC ≥ 7.5) the f-on-GPU path. Not yet exposed in Python: forward-mode sensitivities, time-dependent RHS, and the CppAD auto-diff path.

Usage

User ODE formulation

• CppAD / automatic differentiation: Pass an ADFun<double> (tape) describing the RHS — sensitivities/Jacobians can be handled via AD/Internal Numerical Differentiation.
• Analytical derivatives: Use the user-ODE interfaces in include/yabi/user_ode.hpp and the manual sensitivity examples (examples/manual_sensitivity_demo.cpp, examples/manual_sensitivity_pure.cpp) to provide analytic Jacobians/sensitivity RHS without requiring CppAD.
• Custom sparse solvers: Implement yabi::LinearSolver (include/yabi/linear_solver.hpp) to plug your own sparse direct/iterative solver; set it via integrator.setLinearSolver(...). UMFPACK/KLU are provided via SuiteSparse; dense LAPACK is the default.

#include <yabi/bdf_integrator.hpp>

CppAD::vector<CppAD::AD<double>> x(1);
x[0] = 1.0;                       // y(0)
CppAD::Independent(x);

CppAD::vector<CppAD::AD<double>> f(1);
f[0] = -0.5 * x[0];

CppAD::ADFun<double> rhs;
rhs.Dependent(x, f);
rhs.optimize();

yabi::BDFIntegrator<double> integrator(rhs);

double y0[1]       = {1.0};
double abs_error[1] = {1e-12};
double y[1];

integrator.resetIntegration(y0, 0.1, 1e-10, abs_error,
                            false, false, nullptr, nullptr, nullptr, nullptr, 0);
integrator.setAlgorithmSettings(5, std::numeric_limits<double>::epsilon(), true, 1);
integrator.performIntegration(2.0, 100000);
integrator.getSolutionVector(y);

Dense output: getSolutionAt(t, y) evaluates the continuous extension at any t inside the last accepted step [t_end - h, t_end] (it returns false outside):

double h = integrator.getCurrentStepSize();   // last step size
integrator.getSolutionAt(2.0 - 0.5 * h, y);   // interpolate inside the last step

The seven trailing sensitivity-seed arguments can be replaced by a yabi::IntegrationOptions struct (defaults reproduce the false, false, nullptr, …, 0 tail):

yabi::IntegrationOptions<double> opts;             // taping=false, direct=false,
integrator.resetIntegration(y0, 0.1, 1e-10, abs_error, opts);  // seeds=nullptr, numOfderivatives=0

For sensitivity, fill the fields and pass the same struct to any of the three overloads (autonomous / time / time+parameter):

opts.taping = false;
opts.direct = false;
opts.y_seed_list              = yseed_ptr.data();   // double** (per-direction state seeds)
opts.p_seed_list              = pseed_ptr.data();   // double** (per-direction parameter seeds)
opts.relative_error_derivatives = rel_err.data();   // double*  (per-direction rel. tol)
opts.absolut_error_derivatives   = abs_ptr.data();  // double** (per-direction abs. tol)
opts.numOfderivatives         = ndir;
integrator.resetIntegration(y0, t0, p, max_step, 1e-10, abs_error, opts);

The RHS type is inferred from the dimensions of the CppAD tape:

Tape domain	Tape range	Meaning
n	n	autonomous: f(y)
n + 1	n	time dependent: f(t, y)
n + 1 + p	n	time and parameter dependent: f(t, y, p)

setAlgorithmSettings(max_order, machine_precision, hit_end_time, error_weight_index) controls the maximum BDF order (1–6), the precision used for the final-step test, whether the integrator must land exactly on the end time, and which error-weight norm to use. All four arguments are validated before any state is changed.

max_order is a cap, not a forced order: the controller chooses the order actually used. Use getCurrentOrder() to read it. Likewise max_StepLength caps the step but is not an accuracy knob. Control accuracy with rtol/atol; because the controller changes order as the step cap changes, the error is not monotone in max_StepLength (measured examples in tests/system/README.md).

To force a fixed order (e.g. for per-order convergence measurement), use setFixedOrder(q):

integrator.setFixedOrder(4);        // pin the order to 4 (ramps up 1..4, then stays)
integrator.getCurrentOrder();       // read the order actually used

Step-controller adaptation (stability)

The step controller is the original 2011 BDF algorithm with one small, deliberate addition for stability: on a rejected step it estimates the local truncation error at the next-lower order (from the highest Nordsieck term) and reduces the order when that gives a smaller or already-acceptable error. This order reduction lets the integrator recover from stiff-oscillatory problems — e.g. the Oregonator's relaxation oscillations, where high-order BDF is unstable at any step size. It is a no-op for well-behaved problems (they accept on the first attempt, so no order reduction ever fires), and each reduction is recorded in the sensitivity tape so forward/reverse replays reconstruct the correct order sequence.

Extreme-stiffness behaviour (EMEP)

On the hardest stiff benchmarks (the CWI "Test Set for IVP Solvers" EMEP problem — 66 species, stiffness ratio ~1e14) the plain defaults used to stall with too much steps rejected. The cause was a one-character boundary bug in the step controller.

The minimal-step floor sets h = 4·eps·|t| exactly (eps = machine precision), yet the "accept the step at minimal size" escape tested h < 4·eps·|t|. The strict inequality never matched the floored value, so once a stiff transient drove the step down to the floor, the local-error test could no longer be satisfied and the step was rejected forever. The fix tests h <= 4·eps·|t| (and h > 4·eps·|t| in the still-reducible branch), so at the floor the step is accepted with a warning rather than rejected indefinitely — the standard hmin handling (Hairer & Wanner §III.7).

With that, EMEP matches its CWI reference to ~8e-6 with default Newton settings and no fixed order — only the standard stiff-BDF max order of 5 (BDF6 is not A(α)-stable; see tests/system/45_emep.cpp).

setNewtonLineSearch(true) remains an independent, opt-in feature: a damped Newton (backtracking line search) that halves the Newton correction while the full step would increase the BDF residual — the standard fix for Newton overshoot (Nocedal & Wright §11; Hairer & Wanner §IV.8). It is a mathematical no-op where Newton does not overshoot (the full step is accepted), so it is off by default — the one extra residual evaluation reorders floating-point, so the BASELINE anchors stay byte-identical only with it off. EMEP does not need it; it is kept for general Newton-overshoot robustness.

References: Hairer, Nørsett & Wanner, Solving ODEs I/II; Gear (1971), Numerical Initial Value Problems in ODEs; Byrne & Hindmarsh (1975); Nocedal & Wright, Numerical Optimization §11; Mazzia et al., "Test Set for IVP Solvers".

Relationship to the dissertation (method provenance)

The mathematical formulation is the dissertation by Dominik Skanda (University of Freiburg). That method is correct and remains the baseline: a variable-order, variable-step BDF corrector with a Milne (corrector–predictor) local-error estimator and an error-per-step (EPS) acceptance test — the dissertation's original estimator, which is kept.

Two deliberate, literature-referenced robustness additions sit on top of that baseline. Both are no-ops on well-behaved problems (they fire only on the pathological cases they exist for), so the BASELINE anchors are byte-identical to the dissertation method for the standard benchmarks:

1. Order reduction on a rejected step — [§ "Step-controller adaptation" above]. Source: Gear (1971); Byrne & Hindmarsh (1975).
2. Damped Newton (backtracking line search, opt-in) — [§ "Extreme-stiffness behaviour" above]. Source: Nocedal & Wright, Numerical Optimization §11; Hairer & Wanner §IV.8.

A third change is a bug fix, not an addition: the "accept the step at minimal size" escape that was present in the original code (since the initial commit) tested h < 4·eps·|t| while the minimal-step floor sets h = 4·eps·|t| exactly. Correcting the test to h <= 4·eps·|t| makes the escape fire as it was always intended, which is what lets the extreme-stiffness EMEP problem complete (see above). It is likewise a no-op on well-behaved problems: a healthy step is far above the floor, so the comparison outcome is unchanged.

Sensitivity API

Forward mode

Enable taping in resetIntegration, integrate to the output time, then call forward:

double y_seed[2] = {1.0, 0.0};    // d/dy0
double p_seed[2] = {0.0, 0.0};
double direction[2];

integrator.forward(1, y_seed, p_seed, false);
integrator.getForwardSolutionVector(direction);

For a parameter derivative, set y_seed to zero and the corresponding entry of p_seed to one. Calling forward(..., true) at the same output time resets the forward replay so that a new direction can be requested.

Reverse / adjoint mode

Integrate to the output time with taping enabled, then seed the adjoint with the weights of the scalar output of interest:

double weight[2] = {1.0, 0.0};
double* adj_species[2];
double* adj_parameter[2];

integrator.reverse(1, weight);
integrator.getReverseSolutionVector(adj_species, adj_parameter);

adj_species[i][k] is the k-th order adjoint of state i; adj_parameter[i][k] is the corresponding parameter adjoint.

Continuous adjoint (first order, stable reverse)

The discrete adjoint reverse() above is exact for the discretization but A-stability-bound over long stiff intervals (see T24). continuousAdjoint() integrates the adjoint ODE dλ/dt = -J(y)ᵀλ backward in reversed time with the same BDF integrator, so it stays stable where the discrete replay diverges (e.g. HIRES to steady state):

double weight[2] = {1.0, 0.0};   // d g / d y(T)  for  g = w^T y(T)
double adj_y0[2];
double adj_p[1];                 // one entry per parameter (RHS_type 2)
integrator.continuousAdjoint(weight, adj_y0, adj_p);
// adj_y0[i] = d g / d y0_i,   adj_p[j] = d g / d p_j

It requires only that performIntegration has run (the forward trajectory is recorded automatically) and works with both CppAD (AD Jᵀλ/(∂f/∂p)ᵀλ via Reverse) and the pure UserODERHS backend (analytic Jacobian). The parameter sensitivity ∂g/∂p is accumulated with an augmented backward state (the μ accumulator), and the pure backend supports time- and parameter-dependent RHS via setTimeDerivative/setParameterDerivative.

The second-order continuous adjoint is secondOrderAdjoint — it returns the Hessian-vector product H g · s of g = wᵀ y(t_end) in a seed direction s, stably for long/stiff problems (the stable complement to reverse(2)):

double weights[2] = {1.0, 0.0}, direction[2] = {1.0, 0.0}, hessian[2];
integrator.secondOrderAdjoint(weights, direction, hessian);
// hessian[i] = sum_j d^2 g / d y0_i d y0_j * direction[j]

It needs a Hessian-in-direction hook: CppAD provides it via a Reverse(2) sweep automatically; the pure backend adds setHessianVectorProduct(hv).

The second-order forward is secondOrderForward — it returns the second-order directional derivative d²y/dτ² in a seed direction, stably for long/stiff problems (the stable complement to forward(2)):

double direction[2] = {1.0, 0.0}, second[2];
integrator.secondOrderForward(direction, second);
// second[i] = d^2 y_i / d tau^2  along `direction`

It needs a second-order directional-derivative hook: CppAD provides it via a Forward(2) sweep automatically; the pure backend adds setSecondOrderDirectional(sod).

Higher-order forward sensitivity is forwardOrder(k, seed, result) — it returns d^k y/dτ^k for arbitrary k >= 1, stably for long/stiff problems (the stable complement to forward(k), which is A-stability-bound at order >= 3):

double direction[2] = {1.0, 0.0}, kth[2];
integrator.forwardOrder(3, direction, kth);
// kth[i] = d^3 y_i / d tau^3  along `direction`  (k = 2 coincides with secondOrderForward)

For k >= 2 it needs the k-th-order directional source of the RHS via AbstractRHS::orderKDirectional: CppAD provides it automatically (full-Taylor Forward); the pure backend adds setOrderKDirectional(okd) and forwardOrder fails loudly without it. No finite differences anywhere — the pure hook is exact (analytic), and tools/generate_derivatives.py can emit it for you (see below).

Generating higher-order derivatives (pure mode, optional tool)

tools/generate_derivatives.py is a decoupled, optional SymPy code generator for the exact higher-order sources. It is not a build dependency and the solver never calls it — it emits plain C++ in the UserODERHS hook signatures, which you wire in by hand:

python3 tools/generate_derivatives.py            # print the generated header
python3 tools/generate_derivatives.py --check    # verify the sources numerically

#include "generated_myode.hpp"   // the emitted header
yabi::UserODERHS<double> rhs(n, generated::f, generated::jacobian);
rhs.setOrderKDirectional(generated::orderKDirectional);

Edit the N / Y / F / MAX_ORDER block at the top of the script to describe your ODE, then regenerate. Requires SymPy (pip install --user --break-system-packages sympy; see PEP 668). tests/system/generated_cubic.hpp is a committed example of its output, exercised by tests/system/39_generated_higher_order.cpp.

Simultaneous (variational) sensitivities

Pass seed directions and error controls directly to resetIntegration:

double* y_seeds[ndir];
double* p_seeds[ndir];
double* abs_seeds[ndir];
double rel_seed = 1e-9;
// fill y_seeds[d], p_seeds[d] and abs_seeds[d] ...

integrator.resetIntegration(y0, t0, p, max_step, rel_error, abs_error,
                            false, false, y_seeds, p_seeds, &rel_seed,
                            abs_seeds, ndir);
integrator.performIntegration(t_end, max_steps);
integrator.getSimultaneousDerivativesVector(solution_derivatives);

examples/sensitivity_demo.cpp demonstrates all three modes and checks them against central finite differences.

The forward/simultaneous sensitivity directions are independent, so their Newton solve + update parallelize across threads (T35):

integrator.setSensitivityThreads(4); // 1 = serial (default); >1 = parallel

The direction residual (CppAD Forward(1) tape replay) stays serial — the tape is stateful — and the parallel solve requires a thread-safe linear solver (the default dense LAPACK solve is). tests/component/parallel_sensitivity.cpp and tests/system/48_parallel_sensitivity.cpp verify parallel == serial bit-for-bit.

Higher-order (forward / reverse)

forward and reverse take an order argument and support arbitrary order (1, 2, 3, …). Order 1 uses the legacy single-coefficient replay. Orders 2 and up use a full Taylor-polynomial replay of the frozen BDF algorithm: every intermediate quantity is a truncated Taylor polynomial in the sensitivity parameter τ, the nonlinear RHS is evaluated with CppAD's full-Taylor Forward(q, xq) (xq.size() == Domain()*(q+1)), and the returned value is the order-th directional derivative d^order y / d seed^order (i.e. order! times the order-th Taylor coefficient — the pure derivative, not the raw Taylor coefficient).

// A single call computes order q directly; no forward(1), forward(2), ...
// warm-up sequence is needed for q >= 2.
integrator.forward(q, y_seed, p_seed, true);
integrator.getForwardSolutionVector(dir);    // d^q y / d seed^q

integrator.reverse(q, weights);              // q-th-order adjoint
integrator.getReverseSolutionVector(adj_sp, adj_p);   // adj_sp[i][k], k = 0..q-1

reverse(q) returns the adjoint coefficients of orders 1..q: adj_sp[i][k] is the k-th-order directional sensitivity of sum_j weights[j] * y_j with respect to initial state i (and likewise for adj_p). It is aligned to the same pure-derivative convention as forward (each order k is rescaled by k! to undo CppAD's 1/k! Taylor weighting). The reverse is the gradient-level dual of the forward solution-level derivative: for a scalar ODE with unit seed and unit weight, forward(q) and reverse(q)[q-1] coincide.

Limitations:

• forward(order >= 2) requires an RHS backend with full-Taylor support (the CppAD adapter). A UserODERHS-style backend exposes only single-coefficient Forward, so order ≥ 2 fails loudly with YABI_USE_CPPAD=OFF.
• direct sensitivity mode has no full-Taylor state of the implicit solve and stays limited to order 1; higher orders fail loudly.
• reverse(q) must follow the forward(q) that produced the output state (it replays the forward tangent direction), as in the examples above.

The tests/system/28_second_order_sensitivity.cpp driver verifies orders 1–4 against the exact solution of y' = y², and tests/system/33_higher_order_sensitivity.cpp covers a coupled nonlinear system, a non-autonomous system, and a parameter-dependent system.

Choosing a sensitivity method

YABI exposes four pure sensitivity methods; no automatic switching or method mixing happens — the user selects exactly one. They split along two axes: forward vs reverse (what the seed/weight means) and discrete vs continuous (how the sensitivity is propagated).

Method	Direction	Propagation	Stable on long/stiff?
forward(q)	forward (state)	discrete IND (replay of frozen step Jacobians)	order 1 yes; order ≥ 3 diverges
reverse(q)	reverse (adjoint)	discrete IND	order 1 yes; order ≥ 3 diverges
continuousAdjoint / secondOrderAdjoint	reverse	continuous ODE λ' = -Jᵀλ (re-integrated)	yes, orders 1–2
simultaneous seeds / secondOrderForward / forwardOrder(k)	forward	continuous ODE S' = J·S (+ higher-order source)	yes, arbitrary order

The discrete (IND) methods are exact for the discretization and are the fastest and most accurate choice on short intervals — but they replay frozen BDF linearizations, so at order 3–6 over long, stiff intervals the replay is A-stability-bound and can diverge to huge finite values (T24). The continuous methods reformulate the sensitivity as another ODE and hand it to the same stable BDF, inheriting its stability for free; they cost a second (or third) full integration.

Decision table

Your situation	Use
General, short interval, need accuracy + speed	forward(1) / reverse(1) (IND)
Long/stiff interval, forward sensitivity	simultaneous seeds (1st) or secondOrderForward/forwardOrder(k) (2nd+)
Long/stiff interval, reverse sensitivity	continuousAdjoint (1st) or secondOrderAdjoint (2nd)
Higher-order (≥ 3) derivative on a short interval	forward(q) / reverse(q) (full-Taylor replay)
Higher-order (≥ 3) derivative, long/stiff, forward	forwardOrder(k) (continuous, stable)

Backend-agnostic contract — every method works identically under CppAD (YABI_USE_CPPAD=ON) and pure mode (YABI_USE_CPPAD=OFF):

• 1st-order continuous/simultaneous: CppAD provides J (and ∂f/∂p) via AD; pure UserODERHS supplies an analytic Jacobian (plus setTimeDerivative/setParameterDerivative for non-autonomous RHS).
• 2nd-order continuous/simultaneous: CppAD provides H·v via a Reverse(2) or Forward(2) sweep; pure mode supplies the analytic setHessianVectorProduct(hv) (adjoint) or setSecondOrderDirectional(sod) (forward) hooks.
• Higher-order (≥ 3) forward: CppAD provides the k-th-order source via the full-Taylor Forward automatically; pure mode supplies the analytic setOrderKDirectional(okd) hook (or generates it — see above).
• No finite differences are used anywhere. The deliberate gaps are: higher-order continuous adjoint (≥ 3) is not provided (nothing consumes a 3rd-order adjoint), and the discrete higher-order replay (forward(q)/reverse(q), q ≥ 2) requires the CppAD full-Taylor backend and fails loudly in pure mode.

Rejected alternatives (recorded so they are not re-introduced): no auto-patching of IND (a diverged replay is reported, not silently corrected), no on-the-fly conservative-mode switching, and no finite-difference fallback.

Tests

The tests are split into three ctest tiers below tests/:

• tests/unit/ — one focused executable per component (IntegratorConfig, ReplayTape, NordsieckState, StepController, NewtonSolver, SensitivityEngine, AdjointEngine, DenseLUSolver/CSC reference path, SparsityPattern, CscMatrix, CppAD RHS sparsity/J·v/batch J·V, UserODERHS analytic sparsity/J·v/batch J·V, the Krylov solver, the pluggable-preconditioner/forcing contracts, and the sensitivity-linearization seams). These are fast, dependency-light contract tests.
• tests/component/ — cross-component integration: native CSC assembly and the CppAD sparse Jacobian path (native_sparse), matrix-free Newton–Krylov against the direct HIRES trajectory (matrix_free_newton), the same matrix-free path through the pure UserODERHS backend (matrix_free_pure, both configurations), a fully J-free solve with no sparsity pattern and no CSC assembly (matrix_free_jfree, both configurations), plus symbolic reuse with UMFPACK/KLU (suitesparse_sparse, registered only when SuiteSparse is available).
• tests/system/ — the mathematical benchmark suite: linear problems, Robertson, Van der Pol, Brusselator, HIRES, Oregonator, Prothero–Robinson, damped pendulum, E5, large sparse/network systems, forward/reverse/ simultaneous sensitivities, higher-order (forward/reverse, orders 1–4), dense output, fixed order, and failure/tolerance regressions.

Every test checks its result against an independent reference solution. Run a single tier with ctest --test-dir build -L unit, -L component, or -L system; ctest also registers the three self-checking examples under the example label (example_basic, example_nonautonomous, example_sensitivity_demo).

See tests/system/README.md for the measured convergence slopes; the per-order convergence test for orders 2–6 is now expressible via setFixedOrder + getCurrentOrder (tests/system/32_fixed_order.cpp).

License

Apache License 2.0. See LICENSE and NOTICE.

API Quick Start

Basic integration (autonomous)

#include <yabi/bdf_integrator.hpp>
#include <cppad/cppad.hpp>

int main() {
    CppAD::ADFun<double> rhs;
    CppAD::vector<CppAD::AD<double>> x(1);
    x[0] = 1.0;
    CppAD::Independent(x);
    CppAD::vector<CppAD::AD<double>> f(1);
    f[0] = -x[0];
    rhs.Dependent(x, f);
    rhs.optimize();

    yabi::BDFIntegrator<double> integrator(rhs);
    double y0[1] = {1.0}, abs_err[1] = {1e-10};
    integrator.resetIntegration(y0, 2.0, 1e-10, abs_err, false, false, NULL,NULL,NULL,NULL,0);
    integrator.setAlgorithmSettings(5, 1e-16, true, 1);
    integrator.performIntegration(2.0, 100000);
    double y[1]; integrator.getSolutionVector(y);
    return 0;
}

Forward/Reverse/Simultaneous sensitivities

See examples/sensitivity_demo.cpp for full usage of forward(), reverse(), and simultaneous seeds via resetIntegration with seed arrays.

Pluggable linear solvers (sparse)

#include <yabi/bdf_integrator.hpp>
#include <yabi/detail/umfpack_linear_solver.ipp>

yabi::BDFIntegrator<double> integrator(rhs);
integrator.setLinearSolver(new yabi::UmfpackLinearSolver<double>()); // takes ownership

Defaults to dense LAPACK. Supports UMFPACK/KLU via SuiteSparse. See examples/heat1d_umfpack_literature.cpp.

GPU linear solver (experimental, T36/T36b/T36c)

An opt-in cuSOLVER dense-LU solver offloads the Newton factorize/solve to a GPU through the same setLinearSolver seam:

#include <yabi/bdf_integrator.hpp>
#include <yabi/cuda_linear_solver.hpp>

integrator.setLinearSolver(new yabi::CudaLinearSolver<double>());

Build it with -DYABI_USE_CUDA=ON (default OFF — the regular build stays CUDA-free, with zero CUDA headers/libraries). It uses only the host-side cuBLAS/cuSOLVER C API (no nvcc), and links -lcusolver -lcudart.

Each instance owns its own CUDA stream + cuSOLVER handle (cusolverDnSetStream) plus pinned host staging buffers (cudaMallocHost) with stream-ordered cudaMemcpyAsync copies, so transfer can overlap compute and multiple instances run in parallel without serializing on the default stream (tests/component/cuda_parallel.cpp runs 4 concurrent integrators, each matching the CPU reference). tests/component/cuda_solver.cpp cross-checks the GPU state solve against dense LAPACK on HIRES to ~3e-16.

Scope (honest): this offloads the linear solve only — the Jacobian J is still assembled on the host. The sparse paths are the cuDSS adapter (T36d) and the GPU Krylov solver (T36e, below). Throughput on a small card (GTX 1050) is modest: ~1.2–1.3× for 4 concurrent instances at n ≈ 1000–2000, slower at n ≈ 400 (contention). The cuBLAS/cuSOLVER API is stable across CUDA 12.6 → 13, so the solver is source-identical on both (only -arch/driver change). Requires CUDA (built against 12.6) + an NVIDIA driver.

GPU sparse iterative Krylov solver (own implementation, T36e)

For large sparse Newton systems (n ≫ 100) where the dense cuSOLVER factorization is too costly, CudaKrylovSolver runs a BiCGStab Krylov solve on the GPU — cuBLAS for the vector ops and cuSPARSE SpMV for P·v — behind the same setLinearSolver seam. It is our own implementation over CUDA primitives (no external library):

#include <yabi/cuda_krylov_solver.hpp>
auto *solver = new yabi::CudaKrylovSolver<double>();
solver->setTolerance(1e-10, 0.0, 200, 30);
solver->setForcingTerm(yabi::eisenstatWalkerForcing<double>()); // inexact Newton
integrator.setLinearSolver(solver);

Built with -DYABI_USE_CUDA=ON (links -lcusparse -lcublas -lcudart, no nvcc; the default build stays CUDA-free). tests/component/cuda_krylov_solver.cpp cross-checks it against dense LAPACK on HIRES to ~3e-16; tests/system/cuda_krylov_heat.cpp solves the method-of-lines heat equation (n=500) to the analytic solution.

External sparse solvers (open architecture, T36d)

The solver architecture is open: setLinearSolver takes any LinearSolver, and include/yabi/external_sparse_solver.hpp provides ExternalSparseSolver — the reusable, solver-agnostic half (CSC→CSR handoff + factorize/solve plumbing). A concrete adapter implements only buildCsc (build the solver from CSR) and solveOne (solve one RHS), so adding a solver is a ~30-line adapter that never touches the integrator core. We integrate sophisticated open-source solvers, we do not reimplement them.

The first adapter is AMGCL (MIT, header-only AMG + Krylov):

#include <yabi/bdf_integrator.hpp>
#include <yabi/amgcl_linear_solver.hpp>

integrator.setLinearSolver(new yabi::AmgclLinearSolver<double>()); // BiCGStab + ILU(k)

Build with -DYABI_USE_AMGCL=ON (default OFF — dependency-free). Provision the headers once (no sudo): AMGCL → ~/.local/opt/amgcl, Boost → ~/.local/opt/boost (download the AMGCL tarball and libboost1.83-dev via apt download, extract). tests/component/amgcl_solver.cpp cross-checks it vs dense LAPACK on HIRES to 4.4e-16. The BDF Newton matrix is non-symmetric, so the default is BiCGStab+ILU; AMG+CG is the SPD-only alternative.

The second adapter is Ginkgo (BSD, many-core OpenMP/CUDA/HIP):

#include <yabi/bdf_integrator.hpp>
#include <yabi/ginkgo_linear_solver.hpp>

integrator.setLinearSolver(new yabi::GinkgoLinearSolver<double>()); // BiCGStab (OpenMP)

Build with -DYABI_USE_GINKGO=ON. Provision Ginkgo from its official source (cmake + install to ~/.local/opt/ginkgo, OpenMP backend). It uses the OpenMP (CPU) backend, so it works without a GPU; the CUDA backend is the same adapter with a CudaExecutor (and a supported GPU arch). tests/component/ginkgo_solver.cpp cross-checks it to 1.1e-16.

The third adapter is cuDSS (NVIDIA's sparse direct solver, via the nvidia-cudss-cu12 wheel):

#include <yabi/bdf_integrator.hpp>
#include <yabi/cudss_linear_solver.hpp>

integrator.setLinearSolver(new yabi::CudssLinearSolver<double>()); // general sparse LU

Build with -DYABI_USE_CUDSS=ON. Provision it with pip download nvidia-cudss-cu12 + extract to ~/.local/opt/cudss (never vendored). cuDSS has real ANALYSIS/FACTORIZATION/SOLVE phases, so it matches the BDF factorize-once/solve-many pattern; tests/component/cudss_solver.cpp cross-checks it to 1.7e-18. All three adapters reuse the same ExternalSparseSolver base — the libraries are never vendored, only the thin adapter + CMake detection live in this repo.

Matrix-free Newton–Krylov (state solve)

#include <yabi/bdf_integrator.hpp>
#include <yabi/matrix_free_solver.hpp>

auto *solver = new yabi::GmresLinearSolver<double>();
solver->setPreconditioner(yabi::KrylovPreconditioner::Jacobi); // default
solver->setTolerance(1e-12, 0.0, 200, 30);                     // rtol, atol, max it, restart
integrator.setLinearSolver(solver); // takes ownership; auto-enables matrix-free assembly

FgmresLinearSolver and BicgstabLinearSolver are drop-in alternatives. The Krylov iteration drives the state Newton solve through the RHS Jacobian–vector product (AbstractRHS::jacobianVectorProduct, one CppAD Forward(1) sweep for AD tapes) and never sees a stored dense J. The scaled CSC P is still assembled for the Newton scaling and the preconditioner (Jacobi or the exact Direct solve). The path requires an RHS sparsity pattern; without one the historical dense solve is used and a single warning is logged instead of silently downgrading. operatorApplications() reports how many J·v products ran, so a test can prove the Krylov path was taken.

Without CppAD, supply the same information on UserODERHS explicitly:

#include <yabi/user_ode_rhs.hpp>

yabi::UserODERHS<double> rhs(n, f, jac);
yabi::SparsityPattern pattern;          // column-oriented df/dy structure
pattern.rows = pattern.cols = n;
pattern.column_rows = {{0, 1}, {0, 1}}; // per-column row indices
rhs.setSparsityPattern(pattern);        // validated; throws if malformed
rhs.setJacobianVectorProduct([](const double *y, const double *v, double *Jv) {
    // analytic J(y)*v; J is never formed
});

Both hooks are additive: without them UserODERHS keeps the historical dense behaviour. This activates the same CSC sparse path used by CppADRHS, in YABI_USE_CPPAD=OFF builds too.

Pluggable preconditioners, J-free solves, and inexact-Newton forcing

The built-in preconditioners remain None, Jacobi, and Direct, but a user object now takes precedence. Implement yabi::Preconditioner<Scalar> (or wrap callbacks with yabi::FunctionPreconditioner) and install it on the solver:

#include <yabi/preconditioner.hpp>

class MyAmg : public yabi::Preconditioner<double> {
  public:
    bool requiresAssembledMatrix() const override { return true; }
    bool update(const yabi::CscMatrix<double> *P,
                const yabi::MatrixFreeOperator<double> *op,
                const std::vector<double> &diagonal, unsigned int n) override {
        // Build/reuse an ILU/AMG/block preconditioner from P (or op).
        return true;
    }
    bool apply(const double *r, double *z) override {
        // z = M^{-1} r.
        return true;
    }
};

MyAmg my_amg;
solver->setUserPreconditioner(&my_amg); // borrowed; caller owns it

This is the adapter boundary for external packages (AMGCL, hypre, Trilinos, block/physics-based solvers); YABI does not vendor ILU/AMG. A preconditioner whose requiresAssembledMatrix() returns false enables a genuinely J-free state solve: no sparsity pattern is requested and no CSC J/P is assembled. update() then receives a null matrix, the live scaled P·x operator, and an empty diagonal, and D1 is estimated with one J·v product in direction 2^D2 (documented as a signed row-sum approximation). A failed Krylov solve reports failure rather than silently falling back, because no assembled matrix exists.

An inexact-Newton forcing term can adjust the relative linear tolerance from the nonlinear residual:

solver->setForcingTerm(yabi::eisenstatWalkerForcing<double>());
// or any callable taking NewtonLinearizationContext<double>

Without a callback the fixed setTolerance() value is used unchanged; lastRelativeTolerance() reports the value actually used.

Batched and matrix-free sensitivity interfaces

AbstractRHS::jacobianMatrixProduct(x, directions, columns, result) computes J·V for several column-major directions, with a correct per-column default, one-Jacobian evaluation for UserODERHS, and one multi-direction CppAD Forward(1) sweep. MatrixFreeOperator::applyBlock is the batched operator hook, and yabi::SensitivityLinearization (with yabi::MatrixFreeSensitivityLinearization as an adapter) names the contract a future sparse/matrix-free sensitivity solve consumes. These are interfaces only: the current forward/reverse/simultaneous replay still uses explicit J/H and its numbers are unchanged.

Sensitivities are deliberately unchanged: forward/reverse/simultaneous sensitivity still uses the explicit J/H blocks and is not matrix-free. tests/component/matrix_free_newton.cpp verifies all three Krylov methods on HIRES against the direct trajectory and checks that the sensitivity result is unchanged; tests/component/matrix_free_pure.cpp does the same through UserODERHS and runs in both build configurations; tests/component/matrix_free_jfree.cpp exercises the no-pattern J-free path.

Selectable Jacobian backends (T33)

How the Newton-matrix Jacobian is formed for the plain state solve is chosen explicitly through setJacobianBackend:

#include <yabi/jacobian_backend.hpp>

integrator.setJacobianBackend(yabi::JacobianBackend::Sparse);  // or Dense / MatrixFree / Auto

• Dense — J/P via Forward(1) directional columns, dense LAPACK. Always available; the default solver path.
• Sparse — J/P in CSC via sparseJacobian(), factorized by a CSC-capable solver. Requires a RHS Jacobian sparsity pattern + supportsCsc() solver.
• MatrixFree — J·v via jacobianVectorProduct(), solved by a Krylov solver (sparse or J-free preconditioner). Requires a supportsMatrixFree() solver, plus a pattern if its preconditioner needs an assembled matrix.
• Auto (default) — the historical capability-driven selection: dense, unless you install a UMFPACK/KLU solver (→ sparse) or a Krylov solver (→ matrix-free).

The Jacobian source is orthogonal and is chosen by which AbstractRHS you construct: CppADRHS (AD) or UserODERHS (analytic callbacks) both feed every representation above, so "user-analytic" is not a separate backend value.

An explicit Dense/Sparse/MatrixFree choice is honoured strictly: if the RHS + installed solver cannot satisfy it, the solve fails loudly rather than silently downgrading to dense. The sparse/matrix-free paths cover the plain state Newton solve only — sensitivity solves (simultaneous seeds, or taped forward/reverse replay) always use the explicit dense J/H path, so nothing is silently lost. getJacobianBackend() reports the selection.

Pure/manual ODE interface (no CppAD dependency in user code)

Users provide ODE as std::function<void(double,const double*,int,double*)> via include/yabi/user_ode.hpp. See tests/system/manual_sensitivity_pure.cpp.

Tests

ctest --test-dir build --output-on-failure   # all tiers
ctest --test-dir build -L unit               # isolated component contracts
ctest --test-dir build -L component          # native sparse, matrix-free, SuiteSparse
ctest --test-dir build -L system             # mathematical benchmarks
./run_tests.sh                                # verbose full run

Detailed usage

1) Providing the ODE (user function)

• Via CppAD tape (AD): Build CppAD::ADFun<double> from RHS f(x) (or f(t,x)) — see Usage example. CppAD optional.
• Via analytical derivatives / manual interfaces: Define RHS and provide analytic Jacobian/derivative routines through include/yabi/user_ode.hpp. See examples/manual_sensitivity_pure.cpp and examples/manual_sensitivity_demo.cpp (no CppAD tape required in pure form).

2) Analytical derivatives

• Supply analytic Jacobian df/dy and/or derivative RHS for sensitivity equations as appropriate for your formulation.
• Use the user-ODE interface; see include/yabi/user_ode.hpp for the API and examples/manual_sensitivity_pure.cpp for a complete example.

3) Pluggable linear solvers (dense + sparse)

Interface: include/yabi/linear_solver.hpp defines the abstract solver used by the corrector/linear solves. Implement initialize(), factorize() and solve(); a solver that consumes native CSC storage additionally overrides supportsCsc(), factorizeCsc() and solveCsc().

Built-in:

• DenseLUSolver<double> (default, LAPACK dgetrf/dgetrs)
• UmfpackLinearSolver<double> (SuiteSparse UMFPACK)
• KluLinearSolver<double> (SuiteSparse KLU)
• GmresLinearSolver<double> / FgmresLinearSolver<double> / BicgstabLinearSolver<double> (matrix-free state Newton solve via J·v; see above)

Connect your own solver (complete example):

#include <yabi/linear_solver.hpp>
#include <yabi/bdf_integrator.hpp>

class MySparseSolver : public yabi::LinearSolver<double> {
public:
  bool initialize(unsigned int n) override { /* size workspace */ return true; }

  // Dense interface: P is column-major n x n.
  bool factorize(double *P, int *ipiv, long *scalingD1, long *scalingD2,
                 unsigned int n) override { /* build symbolic/numeric */ return true; }
  int solve(double *P, int *ipiv, long *scalingD1, long *scalingD2,
            double *b, unsigned int n, unsigned int nrhs) override { /* solve */ return 0; }

  // Optional native CSC path (T31): return true to receive CscMatrix<double>.
  bool supportsCsc() const override { return true; }
  bool factorizeCsc(const yabi::CscMatrix<double>& P, int *ipiv,
                    long *scalingD1, long *scalingD2, unsigned int n) override { return true; }
  int solveCsc(const yabi::CscMatrix<double>& P, int *ipiv,
               long *scalingD1, long *scalingD2, double *b,
               unsigned int n, unsigned int nrhs) override { return 0; }
};

yabi::BDFIntegrator<double> integrator(rhs);
integrator.setLinearSolver(new MySparseSolver()); // integrator takes ownership

See include/yabi/detail/dense_linear_solver.ipp, umfpack_linear_solver.ipp, klu_linear_solver.ipp for reference implementations. Any sparse solver (MUMPS, SuperLU, Pardiso, Cholmod, etc.) can be plugged this way.

Complete examples: examples/heat1d_umfpack_literature.cpp, examples/heat1d_dense_literature.cpp, examples/manual_sensitivity_pure.cpp, examples/manual_sensitivity_demo.cpp.