Prototype: exact Jacobian from Enzyme and the sparse KLU solver for CVODE (#1412, #1409) by matthiaskoenig · Pull Request #1418 · sys-bio/roadrunner

matthiaskoenig · 2026-11-03T06:02:20Z

This draft is a working prototype of option 3 of #1412: CVODE gets an exact Jacobian, which Enzyme computes from roadrunner's generated LLVM IR, and the sparse direct solver KLU (#1409) uses it. It also partly covers #1411: the CVODE counters of a simulation are available from Python.

It is a draft for discussion, not for merging as it is. Everything is opt in, behind two CMake options that are off by default and one load option that is off by default. The questions it answers:

Is it exact? Yes. It matches finite differences on four benchmark models and on all 1417 SBML test suite cases it covers, with MCJit and with LLJit.
Is it fast? One Jacobian is 2 to 50 times faster than CVODE's dense difference quotient. A whole simulation is 1.2 to 5.3 times faster, but only together with KLU. The JIT does not change these ratios.
What does it cost to load? A few seconds for the largest model. Chen_MSB2009 (500 states, 827 reactions) loads in 5.2 s with LLJit, against 1.1 s without the Enzyme Jacobian; with MCJit it is 8.4 s against 3.5 s. The three smaller models load 0.2 to 0.7 s slower.
Does it cost anything when switched off? No. Loading and simulating stay within 3 % of develop.

Everything was measured with both JITs. Measuring with LLJit showed that most of the load time went into the Enzyme pass, not into compiling the derivative code. Differentiating the reaction rates in blocks of 64 reactions cuts the pass to a third; Chen_MSB2009 took 38 s to load in the first MCJit version. The full write-up, with the tables, is prototype/enzyme_jacobian/RESULTS.md.

Background

Why a Jacobian. CVODE's stiff solver (BDF with Newton iterations) needs the Jacobian J = ∂f/∂y of the right-hand side. roadrunner does not provide one, so CVODE estimates it with a dense difference quotient (cvLsDenseDQJac). That costs one evaluation of the model per state, and it has a finite difference error. In the profiles of #1409, cvLsDenseDQJac takes 10 to 15 % of the instructions of a simulation. More importantly, every sparse linear solver in SUNDIALS needs a sparse Jacobian from the user, so without a Jacobian from roadrunner there is no KLU. The Jacobian would also serve getFullJacobian, the steady-state solvers and forward sensitivities (see #1412).

Why Enzyme. #1412 lists the options: symbolic derivatives of the SBML math with libsbml, generated derivative code, or automatic differentiation. Enzyme is an LLVM plugin that differentiates LLVM IR at compile time. roadrunner already turns the model into LLVM IR, so Enzyme can differentiate exactly the code CVODE evaluates:

The generated evalReactionRates and evalRateRuleRates are differentiated as they are.
Assignment rules, function definitions, piecewise, events and conversion factors come along without a second implementation of the math.
The derivative is exact to rounding error.

Enzyme's forward mode (__enzyme_fwddiff) computes a Jacobian-vector product J·v in about the cost of one evaluation of the model. With column coloring, one such pass gives a whole group of columns. Here that is 2 to 7 passes per Jacobian instead of n model evaluations.

The go/no-go cutoff we agreed on for the prototype: it was viable only if Enzyme ran against roadrunner's prebuilt LLVM without rebuilding LLVM. It does. The tagged Enzyme release v0.0.303 works against roadrunner's LLVM 23.1.4 (static, RTTI, no clang) without patches, built as the static library libEnzymeStatic-23.a and run in the new pass manager inside roadrunner's JIT.

How it works

Build. -DRR_USE_ENZYME=ON -DRR_ENZYME_LIBRARY=.../libEnzymeStatic-23.a links Enzyme. -DRR_USE_KLU=ON -DRR_SUITESPARSE_PREFIX=... links KLU from a SuiteSparse install and needs SUNDIALS built with KLU. Both options are off by default.
Load. Config.LLVM_ENZYME_JACOBIAN (load option ENZYME_JACOBIAN) turns the feature on for a model. After code generation, EnzymeJacobianCodeGen::prepare:
- computes a conservative sparsity pattern from the SBML math (JacobianSparsity). It walks the ASTs through assignment rules, function definitions and all branches of piecewise, and handles local parameters, compartments of species and rate rules on concentration species;
- colors the columns (greedy distance-1 coloring);
- emits evalReactionRatesTangents and evalRateRuleRatesTangents, which call __enzyme_fwddiff on the generated functions with a shadow LLVMModelData;
- splits evalReactionRates into blocks of 64 reactions. Each block is a copy that keeps only the stores of its reaction rates, with the code they no longer need removed; Enzyme differentiates each copy, and the copies are deleted after the pass. The split requires that the function writes only the reaction rates and never reads them back, which holds because a reaction id in the math is computed again from its kinetic law. This is checked on the IR, and otherwise the function is differentiated whole;
- runs the Enzyme pass.
The model cache key includes the option.
Fallback. If a model cannot be differentiated, the model loads without the Enzyme Jacobian, gives the reason through getStateVectorJacobianUnavailableReason(), and CVODE keeps the difference quotient. Enzyme reports failures as LLVM error diagnostics, which terminate the process by default. A custom DiagnosticHandler turns them into a fallback instead.
Evaluation. LLVMExecutableModel::getStateVectorJacobian(t, y, values):
- evaluates the model once;
- runs one tangent pass per color;
- decompresses the result into dv/dy and d(rate rules)/dy;
- assembles J = cf · N · dv/dy, plus the rows of the rate rules, directly in CSC order.
CVODE. Two new integrator settings, jacobian ("difference_quotient" by default, or "enzyme") and linear_solver ("dense" by default, or "klu"). CVodeSetJacFn gets one callback, which fills the dense matrix or copies the pattern and values into a SUNSparseMatrix for SUNLinSol_KLU. If a combination cannot be served (a model without the Jacobian, or a build without KLU), CVODE falls back to dense with the difference quotient and logs a warning. This also covers loadState, copies and conserved moiety analysis. The one error is jacobian = "enzyme" on a model loaded without the option.

Generated helpers that Enzyme cannot see into (rr_min, rr_max, rr_logd, rr_rootd, arccot, rr_arccot_negzero) get IR bodies in the module.

API

import roadrunner
roadrunner.Config.setValue(roadrunner.Config.LLVM_ENZYME_JACOBIAN, True)
r = roadrunner.RoadRunner("model.xml")
r.integrator.setValue("jacobian", "enzyme")
r.integrator.setValue("linear_solver", "klu")       # needs RR_USE_KLU
r.simulate(0, 100, 101)
r.integrator.getSolverStatistics()
# {'steps': 587, 'rhs_evals': 620, 'rhs_evals_jac': 0, 'jac_evals': 10, 'nonlin_iters': 617,
#  'nonlin_conv_fails': 0, 'err_test_fails': 4, 'lin_setups': 52, 'sparse_linear_solver': 1}
r.getStateVectorJacobian()          # dense matrix over the state vector
r.getStateVectorJacobianPattern()   # its sparsity pattern

C++: ExecutableModel::hasStateVectorJacobian(), getStateVectorJacobianUnavailableReason(), getStateVectorJacobianPattern() (rr::CscPattern) and getStateVectorJacobian(time, y, values), plus CVODEIntegrator::getSolverStatistics(). The state vector is the one CVODE integrates: the rate rule values, then the amounts of the independent floating species.

getSolverStatistics() reports the counters of #1411 under shorter keys and without the totals over reinitializations that #1411 asks for. It was needed to measure this prototype, and it is not a proposal for the final names.

JIT backends

roadrunner has two JITs, and I measured three configurations: mcjit; lljit-none, LLJit as it is configured today; and lljit-aggressive, the codegen level #1215 apparently meant. Three findings about the JITs shaped the measurements:

MCJit is still the default. Change the default backend to LLJit. #1215 changed the map of the configuration file names instead of the default value (MCJit is still the default JIT: #1215 changed the configuration file key map instead of the default, which also breaks three keys of the configuration file #1419).
LLJit runs no IR optimization passes. Only MCJit runs the OPTIMIZE_* passes, so with LLJit the codegen level is the only optimization of the generated code, Enzyme's included.
LLJit compiles a model only once per process. It reuses the object regardless of LOADSBMLOPTIONS_RECOMPILE and of the optimization level (LLJit reuses the compiled object of a model in a process despite LOADSBMLOPTIONS_RECOMPILE and a changed optimization level #1420). The benchmark therefore loads every model in a fresh process and runs every backend in a process of its own.

Results

Four models:

hctz_body and glimepiride_body: whole-body PBPK models, the flat SBML of hctz-model and glimepiride-model.
Smith_BMCSystBiol2013 and Chen_MSB2009: from the PEtab benchmark collection, at their nominal parameters.

Settings: rtol = atol = 1e-10, simulate(0, end, 201). Each backend runs in its own process, and the backends take turns for three rounds. Every time is the median over the rounds of 10 simulations per configuration (the configurations take turns) and 3 loads, each in a fresh process.

model	states	backend	load (s)	load with Enzyme (s)	Jacobian, DQ (µs)	Jacobian, Enzyme (µs)	speedup	simulate, dense + DQ (ms)	dense + Enzyme	KLU + Enzyme
hctz_body	28	mcjit	0.28	0.56	11	5.86	1.9x	2.21	1.03x	1.16x
		lljit-none	0.12	0.30	25.4	12.3	2.1x	2.08	1.07x	1.20x
		lljit-aggressive	0.19	0.39	11.7	5.34	2.2x	2.08	1.07x	1.23x
glimepiride_body	51	mcjit	0.44	0.85	12	1.42	8.5x	2.49	0.95x	1.31x
		lljit-none	0.15	0.54	18.5	2.34	7.9x	3.2	1.06x	1.28x
		lljit-aggressive	0.29	0.70	18.6	1.96	9.5x	1.58	1.06x	1.38x
Smith_BMCSystBiol2013	104	mcjit	0.94	1.68	72.2	4.4	16.4x	19.7	1.16x	1.73x
		lljit-none	0.32	0.87	74.8	4.82	15.5x	22	1.16x	1.74x
		lljit-aggressive	0.64	1.30	69.5	6.13	11.3x	19.1	1.11x	1.70x
Chen_MSB2009	500	mcjit	3.47	8.40	2570	53.1	48.4x	1930	1.05x	4.66x
		lljit-none	1.08	5.23	2826	64.2	44.0x	1820	0.98x	4.80x
		lljit-aggressive	1.70	7.17	2643	51.8	51.0x	2144	1.07x	5.29x

The machine was shared with other jobs (load average 12 to 22 on 8 performance and 4 efficiency cores), and a process on an efficiency core runs 2 to 3 times slower. Single times therefore scatter, for example the lljit-none Jacobian of hctz_body. The ratios within a backend are stable within about 10 %.

What the numbers show:

The sparse linear algebra brings the speedup, and the Jacobian makes it possible. With the dense solver, the exact Jacobian gains at most 16 %. CVODE reuses its Jacobian and evaluates it only 6 to 27 times per simulation, and on these models evaluating the model is cheap compared with the dense LU. KLU factorizes the sparse I - γJ, and the gain grows with the size of the model, as the profiles in Sparse linear solver (KLU) and a sparse Jacobian for CVODE #1409 predicted.
On Chen_MSB2009, KLU also takes fewer steps. KLU + Enzyme takes 587 steps instead of 1349, with 4 error test failures instead of 143. Dense + Enzyme takes about as many steps as dense + DQ, so the difference comes from the linear solver. My hypothesis is that dense LU loses accuracy on this badly scaled matrix (entries up to 1.2e9) and KLU does not; I have not verified it.
The Enzyme pass, not the JIT, was the main load cost. LLJit at none compiles the derivative code of Chen_MSB2009 in about 1 s, while the Enzyme pass on its 827 reactions took 8.6 s. Sampling with gdb put the time into Enzyme's activity and type analysis, which grow faster than linearly with the size of the function (0.14 s for 31 reactions, 1.1 s for 367, 8.6 s for 827). In blocks of 64 reactions the pass takes 2.6 s. MCJit benefits too: its compilation of the smaller derivative functions drops from 10.7 s to 4.9 s. blocks.csv has the sweep over block sizes. Enzyme's command line options did not help.
The JIT matters for the derivative code, not for the ratios. Without IR optimization, LLJit at none makes the derivative code slower (Chen_MSB2009: 64 µs per Jacobian against 52 to 53 µs), which the simulations hardly notice. Enzyme's vector mode (several colors per pass) barely speeds up the Jacobian but multiplies the compile time, so the default width is 1.

Validation

Jacobian. Compared against central differences with a ladder of step sizes, within 1e-6 relative, and one-sided differences within 1e-4 at the kinks of min/max, with MCJit and with LLJit at none and aggressive:
- Benchmark models, at three states each: no mismatches, and no finite difference nonzero outside the pattern. The largest relative error is 1e-9.
- SBML test suite, 1823 semantic cases: 1417 pass, none fail, with every backend. 164 fall back with a reason: 104 have a stoichiometry set by a rule, 46 have species conversion factors, and 14 call functions without a derivative (rateOf through csr_matrix_get_nz, factorial, sec). 242 are skipped: roadrunner cannot simulate 208 of them (algebraic rules, delays, fast reactions), and 34 have no timecourse.
- Rounding noise. In one case (01208, LLJit at aggressive), an entry that is zero by cancellation is 2e-16 in the Enzyme Jacobian and up to 9e-11 in the finite differences, because fused multiply-adds keep the rate from cancelling exactly. The validation accepts such an entry when the finite differences of different steps disagree with each other by more than they differ from the Enzyme value, and lists it.
- Blocks. With one reaction per block, 528 test suite cases are split into 2 to 52 blocks, and all of them still pass.
Simulations. All three configurations stay within the integration error of a reference at rtol = 1e-12, atol = 1e-14.
Tests. There are 29 new C++ tests (EnzymeJacobianTests, JacobianSparsityTests, CvodeEnzymeJacobianTests); the 21 Enzyme and CVODE tests run with both JITs, 50 test runs in all. They cover the gate on a real model, the full Jacobian against finite differences, the helper derivatives, every fallback reason, setValue between Jacobians, the cache key, degenerate models, the independence of the Jacobian from the block size, and the CVODE settings in every order, across loads, copies, conserved moiety analysis, events and reset. The tests on the benchmark models skip when RR_ENZYME_MODELS does not point to the model files.
Bugs found and fixed along the way:
- The model cache key ignored the new option.
- The sparsity pattern missed the amount conversion of rate rules on concentration species (16 test suite cases).
- CVODE threw halfway through create() when a setting could not be served, so the next simulation crashed.

Not covered, and open questions for upstream

Load time. A few seconds remain for a model of 800 reactions. Possible next steps:
- copy only the instructions a block needs, instead of copying the whole function and deleting the rest;
- compile the derivative functions lazily, on the first Jacobian CVODE asks for or in the background (LLJit can do this), so that a model never simulated with it does not pay;
- generate the reaction rates in blocks in the code generator itself.
Which way would you prefer?
Dependencies. Enzyme must match roadrunner's LLVM version exactly, so every LLVM update needs an Enzyme release that supports it. KLU needs SuiteSparse, and a SUNDIALS built with it, in libroadrunner-deps. Both would have to go into the dependency builds and the wheels for all platforms. The prototype has only run on Linux (Ubuntu 24.04, gcc 13.3).
Fallbacks that could be removed. These currently fall back to the difference quotient:
- conserved moiety analysis (differentiate with respect to the independent species);
- stoichiometry set by rules;
- species conversion factors;
- a global conversion factor that depends on the state;
- rateOf (csr_matrix_get_nz could be marked inactive);
- factorial, the distrib functions and the reciprocal trigonometric helpers (IR bodies like rr_min/rr_max);
- models restored with saveState/loadState.
Not used yet. getFullJacobian, getReducedJacobian, NLEQ/KINSOL and forward sensitivities could use the same Jacobian.
Not measured. A sparse difference quotient with the same coloring would also feed KLU, without Enzyme: one model evaluation per color, plus the finite difference error. It might be the cheaper route to most of the KLU speedup, if exactness turns out not to matter.
API. jacobian and linear_solver exist in every build, so they appear in getSettings() and in saved states.
Smaller points left for a real implementation. The helper bodies are not removed again when the Enzyme pass fails. There is no null check after looking up the tangent functions. The fallback is decided by comparing reason strings. The list of LLVMModelData fields that the shadow rebases has no static_assert.

Issues found along the way

These are outside the prototype, and I opened them separately:

Python: a model id named like a RoadRunner method (e.g. reset) replaces that method for every instance and removes it for the rest of the process #1415: a model id named like a RoadRunner method (e.g. reset, in 18 test suite cases) replaces that method for every instance and removes it for the rest of the process.
python_tests_RunStochasticTestSuite finds the stochastic test suite only at fixed relative paths #1416: python_tests_RunStochasticTestSuite finds the test suite only at fixed relative paths.
ModelEditingTests and StateSavingTests fail under ctest -j because they share one output directory #1417: about half of ModelEditingTests/StateSavingTests fail under ctest -j because they share one output directory.
MCJit is still the default JIT: #1215 changed the configuration file key map instead of the default, which also breaks three keys of the configuration file #1419: MCJit is still the default JIT, and three keys of the configuration file map to the wrong settings.
LLJit reuses the compiled object of a model in a process despite LOADSBMLOPTIONS_RECOMPILE and a changed optimization level #1420: LLJit reuses the compiled object of a model within a process, despite LOADSBMLOPTIONS_RECOMPILE and a changed optimization level.

Layout and how to reproduce

source/llvm/EnzymeJacobianCodeGen.*, EnzymeShadow.* and JacobianSparsity.* hold the code generation, the shadow model data and the sparsity. The Jacobian is evaluated in LLVMExecutableModel, and CVODE uses it in CVODEIntegrator.
prototype/enzyme_jacobian/ has the rest:
- BUILD.md: building Enzyme, SuiteSparse, SUNDIALS with KLU, and roadrunner;
- validate.py, with its output in validation.md;
- benchmark.py and the C++ timing tool jacobian_timing.cpp, with their output in jit.csv, baseline.csv and width.csv; blocks.py, with its output in blocks.csv. They take the JIT from RR_LLVM_BACKEND and RR_LLJIT_OPT;
- RESULTS.md.
This directory would not go upstream as it is.
The commits follow the steps: build option, gate, sparsity, full Jacobian, CVODE, KLU, validation, benchmark, review fixes, blocks and the measurement with both JITs.

Builds.

With both options on (and BUILD_PYTHON=ON), all new tests pass with both JITs, and so do the other 3315 tests of ctest. Not counted are the tests that already fail on develop (python_tests_RunStochasticTestSuite finds the stochastic test suite only at fixed relative paths #1416, ModelEditingTests and StateSavingTests fail under ctest -j because they share one output directory #1417).
With both options off and the original libroadrunner-deps, the branch builds. The Enzyme pass and KLU are compiled out, and Config.LLVM_ENZYME_JACOBIAN falls back with the reason "roadrunner was built without Enzyme".

…BIAN

…BML math

…obian settings

…ested with MCJit and LLJit Enzyme's activity and type analysis grow faster than linearly with the size of the differentiated function: 8.6 s for the 827 reactions of Chen_MSB2009 in one function. evalReactionRates is therefore differentiated in copies that each keep only the stores of a block of reaction rates (checked on the IR: the function writes only the reaction rates and never reads them), which cuts the Enzyme pass to a third. The copies are deleted after the pass. RR_ENZYME_BLOCK_SIZE, RR_ENZYME_OPTIONS and RR_ENZYME_DUMP_IR are knobs for experiments. Every Enzyme and CVODE test now runs with both JITs, and the load log measures the compilation up to the lookup of the functions, where LLJit compiles.

validate.py, benchmark.py and enzyme_jacobian_timing take the JIT from RR_LLVM_BACKEND and RR_LLJIT_OPT. LLJit keeps the compiled object of a model for the whole process (sys-bio#1420), so every load is measured in a fresh process and every backend in a process of its own (benchmark.py --backends, jit.csv). blocks.py measures the load time against the block size. The validation accepts entries at the rounding noise of the finite differences and reports the error only where the finite differences agree with each other.

matthiaskoenig · 2026-11-03T09:30:49Z

Update: measured with MCJit and LLJit, reaction rates differentiated in blocks

The first version was measured only with MCJit. I have now measured everything with both JITs, which changed the picture of the load time. Two new commits:

1968d6f: Enzyme differentiates the reaction rates in blocks of 64 reactions. Every Enzyme and CVODE test now runs with both JITs.
e4a4507: the scripts take the JIT from the environment, with new results in jit.csv and blocks.csv and a rewritten RESULTS.md.

The description above is updated as well.

The Enzyme pass was the load cost, not the JIT

LLJit compiles the derivative code of Chen_MSB2009 (500 states, 827 reactions) in about 1 s, but the Enzyme pass itself took 8.6 s. Sampling with gdb put the time into Enzyme's activity and type analysis. That analysis grows faster than linearly with the size of the differentiated function:

model	reactions	Enzyme pass
hctz_body	31	0.14 s
Smith_BMCSystBiol2013	367	1.1 s
Chen_MSB2009	827	8.6 s

Enzyme's command line options did not help.

evalReactionRates is now differentiated in blocks:

Each block is a copy of the function that keeps only the stores of its reaction rates. The code those stores no longer need is removed.
Enzyme differentiates each copy, and the copies are deleted after the pass.
The split requires that the function writes only the reaction rates and never reads them back. That holds because a reaction id in the math is computed again from its kinetic law. The code checks this on the IR; otherwise it differentiates the function whole.

The table shows Chen_MSB2009 for a range of block sizes, as medians of 3 loads in fresh processes:

JIT	reactions per block	blocks	Enzyme pass (s)	compilation (s)	load with Enzyme (s)
MCJit	all	1	11.5	10.7	23.2
MCJit	64	13	3.2	4.9	9.4
LLJit (`none`)	all	1	8.6	0.9	10.1
LLJit (`none`)	256	4	4.0	1.0	5.8
LLJit (`none`)	64	13	2.6	1.0	4.7
LLJit (`none`)	16	52	3.9	0.9	8.0

64 reactions per block is now the default. Smaller blocks cost more in copying and per-call overhead than they save. The time of one Jacobian does not depend on the block size.

Benchmark with both JITs

Three configurations:

mcjit;
lljit-none: LLJit as it is configured today;
lljit-aggressive: the codegen level Change the default backend to LLJit. #1215 apparently intended.

Settings and method:

Each JIT runs in a process of its own, and the JITs take turns for three rounds.
The simulation settings are rtol = atol = 1e-10 and simulate(0, end, 201).
Each time is a median over the rounds.

model	JIT	load (s)	load with Enzyme (s)	Jacobian, DQ (µs)	Jacobian, Enzyme (µs)	simulate, dense + DQ (ms)	dense + Enzyme	KLU + Enzyme
hctz_body	mcjit	0.28	0.56	11	5.86	2.21	1.03x	1.16x
	lljit-none	0.12	0.30	25.4	12.3	2.08	1.07x	1.20x
	lljit-aggressive	0.19	0.39	11.7	5.34	2.08	1.07x	1.23x
glimepiride_body	mcjit	0.44	0.85	12	1.42	2.49	0.95x	1.31x
	lljit-none	0.15	0.54	18.5	2.34	3.2	1.06x	1.28x
	lljit-aggressive	0.29	0.70	18.6	1.96	1.58	1.06x	1.38x
Smith_BMCSystBiol2013	mcjit	0.94	1.68	72.2	4.4	19.7	1.16x	1.73x
	lljit-none	0.32	0.87	74.8	4.82	22	1.16x	1.74x
	lljit-aggressive	0.64	1.30	69.5	6.13	19.1	1.11x	1.70x
Chen_MSB2009	mcjit	3.47	8.40	2570	53.1	1930	1.05x	4.66x
	lljit-none	1.08	5.23	2826	64.2	1820	0.98x	4.80x
	lljit-aggressive	1.70	7.17	2643	51.8	2144	1.07x	5.29x

What changed compared with the first version:

Loading Chen_MSB2009 with the Enzyme Jacobian now takes 5.2 s with LLJit and 8.4 s with MCJit. The first MCJit version took 38 s. The three smaller models load 0.2 to 0.7 s slower than without the Enzyme Jacobian.
The speedups do not depend on the JIT. One Jacobian is 2 to 50 times faster than CVODE's dense difference quotient, and a simulation with KLU is 1.2 to 5.3 times faster.
Without IR optimization, LLJit at none makes the derivative code somewhat slower. On Chen_MSB2009 one Jacobian takes 64 µs, against about 52 µs with the other two configurations; the simulations hardly notice.
Noise. Other jobs shared the machine (8 performance and 4 efficiency cores), so single times scatter; the ratios within a JIT are stable within about 10 %.

Validation with both JITs

All three configurations give the same result: the 4 benchmark models and 1417 SBML test suite cases pass, and none fail.

Rounding noise. In one case (01208, LLJit at aggressive), an entry that is zero by cancellation comes out as 2e-16 in the Enzyme Jacobian. The finite differences give up to 9e-11 there, because fused multiply-adds keep the rate from cancelling exactly. The validation now accepts such an entry when the finite differences of different steps disagree with each other by more than they differ from the Enzyme value, and it lists the entry.
Blocks. With one reaction per block, 528 test suite cases are split into 2 to 52 blocks, and all still pass. A C++ test checks that the Jacobian does not depend on the block size.

Two findings about the JITs, opened as issues

MCJit is still the default JIT: #1215 changed the configuration file key map instead of the default, which also breaks three keys of the configuration file #1419: MCJit is still the default. Change the default backend to LLJit. #1215 changed the map of the configuration file names instead of the default value. The same three lines make LLVM_BACKEND, LLJIT_OPTIMIZATION_LEVEL and LLJIT_NUM_THREADS in a configuration file set other settings. The LLJit defaults actually in effect are codegen level none and one thread.
LLJit reuses the compiled object of a model in a process despite LOADSBMLOPTIONS_RECOMPILE and a changed optimization level #1420: LLJit reuses the compiled object of a model for the whole process, despite LOADSBMLOPTIONS_RECOMPILE and a changed optimization level. A benchmark that loads a model twice in one process does not measure the second compilation. benchmark.py therefore loads every model in a fresh process.

matthiaskoenig · 2026-11-03T09:39:41Z

@luciansmith This looks super promising! I see speedups in integration time with Enzyme Jacobian and KLU throughout all models with ~4-5x speedup in integration for larger models (~500 states) with speedup increasing with model size. Our models will have more states 2000 - 50000 so I expect major speedups 10x - 20x here.

When Jitting with LLJit the times are reasonable even for the larger models with a few seconds overhead for calculation of the Jacobian once (which can be cached).

Also this should not be too complicated to implement (see the prototype in this pull request). This can be an alpha/beta feature which can be activated on demand (with default being no Jacobian and dense matrices).

matthiaskoenig added 10 commits November 3, 2026 01:13

Optional Enzyme build behind RR_USE_ENZYME

97b0dff

Enzyme directional derivative of evalReactionRates behind ENZYME_JACO…

57db6e8

…BIAN

Conservative Jacobian sparsity pattern and column coloring from the S…

66e77fa

…BML math

Full Enzyme Jacobian in vector mode with fallbacks

b991f0e

CVODE takes the Enzyme Jacobian and reports its counters

5d404bc

Sparse KLU solver for CVODE with the Enzyme Jacobian

4e2dab1

Validation of the Enzyme Jacobian against finite differences

13ac53e

Benchmark of the Enzyme Jacobian with dense and KLU solvers

5d101d1

CVODE falls back instead of failing when a model cannot serve the Jac…

c4adeff

…obian settings

Prototype notes link the upstream issues found along the way

4179f1d

matthiaskoenig added 2 commits November 3, 2026 10:21

This was referenced Nov 4, 2026

Analytic Jacobian of the model for CVODE, the steady-state solvers and getFullJacobian #1412

Open

Prototype: sparse difference quotient Jacobian for KLU, without Enzyme #1428

Open

Navigation Menu

Repository navigation

Prototype: exact Jacobian from Enzyme and the sparse KLU solver for CVODE (#1412, #1409) - #1418

Prototype: exact Jacobian from Enzyme and the sparse KLU solver for CVODE (#1412, #1409)#1418
matthiaskoenig wants to merge 12 commits into
sys-bio:developfrom
matthiaskoenig:enzyme-jacobian

matthiaskoenig commented Nov 3, 2026 •

edited

Loading

Uh oh!

matthiaskoenig commented Nov 3, 2026

Uh oh!

matthiaskoenig commented Nov 3, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

matthiaskoenig commented Nov 3, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Background

How it works

API

JIT backends

Results

Validation

Not covered, and open questions for upstream

Issues found along the way

Layout and how to reproduce

Uh oh!

matthiaskoenig commented Nov 3, 2026

Update: measured with MCJit and LLJit, reaction rates differentiated in blocks

The Enzyme pass was the load cost, not the JIT

Benchmark with both JITs

Validation with both JITs

Two findings about the JITs, opened as issues

Uh oh!

matthiaskoenig commented Nov 3, 2026

Uh oh!

This branch has not been deployed

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

matthiaskoenig commented Nov 3, 2026 •

edited

Loading