dense: read a row the Newton matrix holds alone as its own equation by wshlavacek · Pull Request #1471 · lanl/bngsim

wshlavacek · 2026-11-03T00:54:18Z

Fixes the x86 path-agreement DIFF on MODEL1812040005|reorder (25.6x the trajectory's peak, run 41491718656). It also fixes BIOMD0000000908's reversed run: silently wrong on the Linux nightly, a failure on macOS.

What was silently wrong on main

model, dense, species reversed	main x86-64	main arm64	exact	this PR
MODEL1812040005 (t 0–100)	P → 2.6e8, 25.6x the peak away	P = −8e-12	P ≡ 0	P ≡ 0 on both
BIOMD0000000908 (t 0–300)	Linux nightly: S grows, the run completes 1.09e134 of the peak away; under Rosetta, `CV_CONV_FAILURE` at t = 145	`CV_CONV_FAILURE` at t = 149.7	S ≡ 0	agrees with the file's order to 1.4e-13 of the peak
P' = P, V' = 1 − V + 1e8·P (new test)	—	V(40) = 0.992	V = 1 − e^−t, P ≡ 0	exact

In each model, the file's own order and KLU keep the 0 exactly. libRoadRunner keeps MODEL1812040005's P at 0 in both orders. It keeps BIOMD0000000908's S at 0 only in the file's order; reversed, its S reaches 6.8e17.

Why

MODEL1812040005's rate for P is P·g(P, V). With P(0) = 0 the solution is unique: P stays at exactly 0. The a·P^z in V's rate does not change that, because it feeds only V. But g > 0 on t ∈ [1.72, 9.53], at up to +9.4 per unit time, so P = 0 is unstable there: any departure grows by e^61.2 ≈ 3.9e26.

While P = 0, P's row of the Newton matrix has nothing off its diagonal, so x_P = b_P / a_PP = 0 exactly. V's rate reads √P, though, so P's column carries a huge entry on V's row.

File order: P's column comes last, its row is never a pivot elsewhere, and x_P is 0.
Reversed order: P's column comes first and partial pivoting takes V's row as its pivot. The back-solve then computes x_P as the remainder of sums that cancel, and their rounding (about 1e-16 by t ≈ 1–2) is left in it. The sign of that residue decided the answer:
- x86-64: it was positive and grew to P's carrying capacity.
- arm64: it was negative, P^0.5 was NaN, and sbml: reversing listOfSpecies makes MODEL1812040005 fail (P below 0, a*P^z NaN) where the order written runs; RoadRunner runs both #1340's band read P as 0, holding it at −8e-12.

KLU's block-triangular form takes such a row as a block of its own, which is why KLU never had this problem. BIOMD0000000908's S = 0 (S' = +0.277·S) is the same case. PR #1427 explained its reorder DIFF as a fragile model's; the dense solve is what broke it.

The fix

Every dense route now uses bngsim's own dense solver (src/lapack_dense_linsol.cpp). That covers CVODE's, the steady-state solves and the step check's scratch solve.
- Where no BLAS is linked, or the BLAS factor is not asked for, it factors with SUNDIALS' built-in GETRF and back-solves with denseGETRS, as SUNLinSol_Dense does.
- SUNDIALS and LAPACK are untouched.
- The GH fix(codegen): walk the derived-parameter DAG instead of flattening it (#99) #132 accessors still read −1 for this non-adaptive mode, so those contracts and SolverStats are unchanged.
At each factorization, before GETRF overwrites the matrix, it notes the rows with no nonzero off the diagonal and their diagonal. The scan is one branchless pass that the compiler vectorizes. A NaN is not a zero.
After each solve, it sets each of those components to b_i / a_ii, the exact answer of the factored system. On the built-in factor, every other component is SUNDIALS' dense LU's to the bit, and so is every component of a matrix with no such row. On the BLAS factor the other components are dgetrf's, as before.
Not covered: species at 0 that feed each other, a chain (Q' = kQ + P) or a cycle (P' = k(P + Q), Q' likewise), are still solved by the pivoted LU. See "Left alone".

Evidence

Tests fail on main and pass here.

python/tests/test_dense_solve_reads_rows_held_alone.py (12 cases): 6 fail on main and all pass here, on arm64 and on x86-64 under Rosetta.
tests/test_lapack_dense_linsol.cpp: both new row-held-alone cases fail on main (8/10). 10/10 pass here with Accelerate. With BNGSIM_ENABLE_LAPACK_DENSE=OFF (from the review), the cases pass or skip by name, with no warnings.

Path agreement on x86-64 (Rosetta build of this branch, with KLU):

MODEL1812040005|reorder is PASS at 1.3e-12 (baseline EXCEPTION, nightly run 41491718656 DIFF 25.6).
BIOMD0000000908|reorder is PASS at 2.1e-13 (baseline DIFF 1.09e134).
BIOMD0000000908|stable stays "contested", which is right: a start moved off 0 still grows.

Corpus sweep, all 1,309 SBML ODE jobs of rr_parity, default route, main against this branch:

File order:
- 1,172 jobs are bit-identical; 76 differ by at most 1e-6 of the peak.
- 8 differ by 1.2e-6 to 3.6e-5. At 1000x tighter tolerances, main and the branch agree to 1e-8 or better on all eight. At the jobs' own tolerances, both sit within 6e-7 to 4e-5 of that answer: 5 slightly farther than main, 3 closer. That is the integrator's own error.
- No job newly fails.
Reversed order (path agreement's reorder transform; 86 models have too few species and reactions for it), counting runs that disagree with their file order by more than 1e-3 of the peak or fail only reversed: 13 on main, 12 here.
- BIOMD0000000908 is fixed. None is new.
- One changed between bad states: MODEL1112110002's reversed run went from finishing at a 100% disagreement to failing with CV_ERR_FAILURE at t = 7. Main's own run at a tenfold tighter tolerance fails the same way. Its GL rate has a pole at GL·(t − 7.2) = −85 that the trajectory crosses near t = 6, so past it no engine has an answer: GL is 1.2e11 in the file's order and 1.5e7 reversed on main. In libRoadRunner at the job's tolerances it is −5.1e6 at t = 7 and reaches −3.45e7 by t = 13.

Performance (CPU time of run(), alternating processes):

Small dense models: BIOMD0000000801 (N = 4, 255k Jacobians) and MODEL0406793751 (N = 55, 28.5k) are within run-to-run noise of main.
Scan cost per factorization (microbenchmark against SUNDIALS' GETRF): 21–29 ns at N = 4; 0.34–0.49 µs at N = 55 (3.5–4%); 12–15 µs at N = 356 (0.4%); 0.12–0.14 ms at N = 1026 (0.2%).
Fill-routed models (BLAS dgetrf, from the review): the first cut's branchy scan cost 7–12%. The branchless scan brings that to +1.0% to +2.3%: egfr_net (N = 356) +2.1%, fceri_ji +2.1%, before_bunching +1.8%, multisite_phos (N = 1026) +2.3%, fceri_fyn +1.0%. Setups, steps and BLAS factorizations are identical to main.

Review

Round 1 (one sub-agent) found:

C++ coverage guard (HIGH): the BLAS half of the new C++ case hid behind a continue on a no-BLAS host, which the coverage guard refuses. It is now its own case with SKIP().
The state-switch continuity test compares the rate-law jump with max|dx/dt| over every species, so one fast unrelated species (for example a large degrading pool) turns a real jump into "continuous". The saltation term is then dropped... #763 test pinned to rounding (HIGH): test_a_fast_loss_does_not_send_the_crossing_back_to_the_old_judgment[1e7] hit The state-switch continuity test compares the rate-law jump with max|dx/dt| over every species, so one fast unrelated species (for example a large degrading pool) turns a real jump into "continuous". The saltation term is then dropped... #763's closer-than-a-step refusal. Whether it fires there turns on rounding: of 17 kf values, main refused 3 on arm64, and the branch refused 4 on arm64 and 4 on x86-64, not the same ones. No build ever gave the old judgment's wrong answer. The test now runs six kf values and asks each for the right columns or that refusal by name, and that at least four run.
Fill-routed slowdown (HIGH): 7–12% from the branchy scan; fixed as above.
Narrower than claimed (CRITICAL, older than this PR): chains and cycles of species at 0 still leak on dense. The claims are narrowed, and the cases are listed below.
Stale comments (LOW): fixed.

Round 2 (a fresh sub-agent) found no CRITICAL or HIGH issues. Its LOW findings were all fixed:

the numbers and wording above;
main's reversed BIOMD0000000908 completes silently wrong on Linux, not with a failure;
libRoadRunner's reversed BIOMD0000000908 is wrong too;
"every other component" holds on the built-in factor only;
stale comments in the header, cvode_simulator.cpp and the C++ test.

It also mapped the #763 test's kf: below 7e7 the refusal flips between builds across a third to a half of kf. The test's six values are 8e7, 1e8, 7e8, 1e9, 1.3e9 and 5e9. They run right on arm64 main, the arm64 branch and the x86 branch. A mutant that brings the old judgment back fails the test: round 2 saw dW = (0, 3, 0) at 1e8, 1e9 and 5e9, and the test accepts no refusal but #763's closer-than-a-step one.

Round 2 also confirmed:

a no-LAPACK build compiles with 0 warnings, and its C++ cases pass or skip by name;
127/127 test_state_switch_sensitivity.py on arm64 and x86;
the fill-routed cost of +1.2% to +2.3%.

Round 3 checked the final text. It found two wording slips in this description, now fixed, and no defect in the code.

Left alone

Chains and cycles of species at 0 on dense (older than this PR, the same on main). These are the review's repros:
- P' = k(P + Q), Q' = k(P + Q), V' = 1 + 1e8(P + Q) − V, with all three starting at 0. In order PQV, max|P| = 9.6e5.
- P' = kP, Q' = kQ + P, V' = 1 + 1e8·Q − V. In orders QPV and PQV, V is off by 6.4e-3 at rtol 1e-8.
- KLU is right in every order. A dense block-triangular solve would cover both shapes, and BIOMD0000000338's XIIa/K/CA (BIOMD0000000338 and 339: KLU and the dense solver take different branches past a step-size collapse (Fg 311 on KLU, 0.73 on dense and libRoadRunner) #1304). That is a design change of its own.
MODEL1112110002 crosses a pole of its GL rate near t = 6 (above). Every engine is past it.
The codegen cache key carries no architecture. On one Mac, an x86-64 process under Rosetta loads an arm64 dylib that a native process cached under the same name, and fails with a dlopen error. It is loud, and unrelated.

Baseline

Path agreement ran on this PR (run 41556140350, x86-64, AMD EPYC 7763). Six rows are copied verbatim from its verdicts.json into differential.json, and the tally is recomputed:

MODEL1812040005|reorder EXCEPTION → PASS and BIOMD0000000908|reorder DIFF → PASS: this fix.
MODEL2001200002|klu_vs_dense and MODEL5952308332|klu_vs_dense EXCEPTION → PASS: ode: KLU's colored difference quotient reads f through the RHS's non-finite retry #1468's, which the last nightly predates. Both pass identically on main.
BIOMD0000000516|steady_vs_run SKIP → PASS and BIOMD0000000510|steady_vs_run PASS → SKIP ("did not converge, residual 6.3e-9"). On these glycolysis models, whether steady_state(method="integration") converges at rtol 1e-6 turns on the run's path, on main as here:
- on arm64 main, 516 converges at rtol 1.05e-6 but not at 0.95e-6 or 1e-6;
- on this branch under x86-64, it fails at 0.9e-6 and 1.05e-6;
- where it fails, the residual stalls at 3e-9 to 8e-9 against a tol of 1e-9.
Several of their species are constant, so their rows of the Newton matrix are held alone. On main the LU leaks rounding into them: Rib5P_g goes from 0.01 to 0.009999999999999993, and Gly_e from 0 to 4.7e-21. That is why this fix moves these rows.

That run's other alert, MODEL1006230077|locale PASS → TIMEOUT, is runner speed. The model spends its whole cap; its base run takes the same 346 steps and 2.0 s on main and here. Every model ran 1.65x as long as on the previous nightly's EPYC 9V45. #1472 has since listed the row in flaky.json. Replaying the run's verdicts against this baseline gives no alert.

Tests

python/tests/test_dense_solve_reads_rows_held_alone.py (new). Cases:
- the two-species unstable zero in both orders, on each Jacobian route and with sensitivities;
- the KLU control;
- MODEL1812040005 and BIOMD0000000908 in both orders, against their file order, with the unstable species at 0 exactly.
- tests/data/BIOMD0000000908.xml (dePillis2013) is new.
tests/test_lapack_dense_linsol.cpp: three new cases.
- A row held alone is read exactly, on the built-in factor and on the BLAS factor. On the built-in factor, every other component is SUNDIALS' to the bit.
- With no row held alone, or with a NaN off the diagonal, the solve is SUNLinSol_Dense's to the bit.
python/tests/test_state_switch_sensitivity.py: the fast-loss case, as above.
python/tests/test_colored_dq_reads_the_clamp.py: test_model9811206584_runs_on_klu_as_on_dense held KLU to dense within 1e-5 of the peak. That failed both Linux legs. At rtol 1e-6, each route is up to 3.1e-5 from the model at rtol 1e-12, where the two routes agree to 3e-12. So the bound measured how alike the two routes' rounding is:
- on arm64 they are 3.1e-6 apart;
- on x86-64 they are 1.8e-5 apart, because this fix's dense run lands nearer the tight run there (9.8e-6 off it in r, against KLU's 2.7e-5).
Both routes are now held to the tight run within 1e-4 of the peak. That passes on arm64 main, on the arm64 branch and under x86-64.
Full suite (-n 6, arm64, merged with main at b1d5f7a): 17,423 passed and 1 failed, test_run_tests_sh, which fails in any scratch worktree and on main there too.

Where row i of the Newton matrix has no nonzero off its diagonal, x_i is fixed by that row alone, b_i / a_ii, and is exactly 0 for a species held at 0. Partial pivoting does not see that: where column i has a larger entry in another row, that row becomes the pivot, and the back-solve computes x_i as the remainder of sums that cancel, leaving their rounding in it. On an exact state that is unstable, the residue grows. MODEL1812040005's P is a rate-rule target whose law is P*g: it starts at 0 and stays there, but g > 0 on t in [1.72, 9.53], where a departure grows by 3.9e26. V's rate reads sqrt(P), so P's column has a huge entry on V's row. In the file's order P's column comes last and P stays 0. Reversed (path agreement's `reorder`), V's row is the pivot. The residue left in P grew to 2.6e8 on x86-64, 25.6x the trajectory's peak. On arm64 it was negative, and #1340's band held P at -8e-12. KLU, whose block-triangular form takes such a row as a block of its own, keeps P at 0 in either order, and so does libRoadRunner. BIOMD0000000908's S = 0 (PR #1427) is the same case: reversed, S grew to 1.09e134 of the peak on the Linux nightly, and the run failed near t = 150 on macOS. The dense solver now notes those rows at each factorization, before the factor overwrites the matrix, in one branchless pass that vectorizes, and the solve reads each of them as b/a in place of what the back-solve left. A NaN off the diagonal is not a zero. On the built-in factor, every other component is SUNDIALS' dense LU's to the bit, and so is every component of a matrix with no such row; on the BLAS factor they are dgetrf's, as before. Species at 0 that feed each other, a chain or a cycle, are not covered: the pivoted LU still solves them on dense. bngsim's own dense solver (lapack_dense_linsol.cpp) now serves every dense route, the steady-state solves and the step check's scratch solve included. Without a BLAS backend, or where the BLAS factor is not asked for, it factors with the built-in GETRF only, as SUNLinSol_Dense does, and the GH #132 accessors still treat it as not the adaptive solver (-1, no reset). test_a_fast_loss_does_not_send_the_crossing_back_to_the_old_judgment pinned kf = 1e7, where #763's closer-than-a-step refusal turns on rounding: of 17 kf, main refused 3 on arm64, this branch 4 on arm64 and 4 on x86-64, not the same ones, and none gave the old judgment's answer. It now runs six kf that ran right on every build and that a mutant bringing the old judgment back gets wrong, and asks each for the right columns or that refusal by name.

…-rows-held-alone

…this PR moves test_model9811206584_runs_on_klu_as_on_dense held KLU to dense within 1e-5 of each species' peak. At the run's rtol of 1e-6 each route is up to 3.1e-5 off the model at rtol 1e-12 (where the two agree to 3e-12), so the bound measured how alike their rounding is: 3.1e-6 apart on arm64, and on x86-64, where this PR's dense solve lands nearer the tight run (9.8e-6 off it in r, KLU 2.7e-5), 1.8e-5 apart, which failed both Linux legs. Both routes are now held to the tight run within 1e-4 of the peak. Path agreement's run 41556140350 on this PR (x86, before #1472 and #961's steady_state merges), copied verbatim: - MODEL1812040005|reorder EXCEPTION -> PASS, BIOMD0000000908|reorder DIFF -> PASS: this PR's fix. - MODEL2001200002|klu_vs_dense, MODEL5952308332|klu_vs_dense EXCEPTION -> PASS: #1468's, which the last nightly predates; both PASS identically on main. - BIOMD0000000516|steady_vs_run SKIP -> PASS and BIOMD0000000510|steady_vs_run PASS -> SKIP (did not converge, residual 6.3e-9). Whether steady_state(method="integration") converges on these glycolysis models at rtol 1e-6 turns on the run's path: on main, 516 converges at rtol 1.05e-6 and not at 0.95e-6 or 1e-6; on this branch under x86-64, not at 0.9e-6 or 1.05e-6. Where it does not, the residual stalls at 3e-9 to 8e-9 against tol 1e-9. The run's other alert, MODEL1006230077|locale PASS -> TIMEOUT, is runner speed, which #1472 now lists in flaky.json. The tally is recomputed; replaying the run's verdicts against this baseline gives no alert.

wshlavacek added 3 commits November 2, 2026 17:54

Merge remote-tracking branch 'origin/main' into fix/dense-solve-reads…

a246808

…-rows-held-alone

wshlavacek merged commit 6af62d2 into main Nov 3, 2026
44 of 45 checks passed

wshlavacek deleted the fix/dense-solve-reads-rows-held-alone branch November 3, 2026 18:35

Montana mentioned this pull request Nov 4, 2026

dense: give exactly 0 to rows at 0 that read only each other (#1485) #1535

Open

Navigation Menu

Repository navigation

dense: read a row the Newton matrix holds alone as its own equation - #1471

dense: read a row the Newton matrix holds alone as its own equation#1471
wshlavacek merged 3 commits into
mainfrom
fix/dense-solve-reads-rows-held-alone

wshlavacek commented Nov 3, 2026 •

edited

Loading

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

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

Uh oh!

What was silently wrong on main

Why

The fix

Evidence

Review

Left alone

Baseline

Tests

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

wshlavacek commented Nov 3, 2026 •

edited

Loading