Skip to content

Keep sparse W sparsity patterns fixed across Jacobian value changes - #4683

Draft
ChrisRackauckas-Claude wants to merge 6 commits into
SciML:masterfrom
ChrisRackauckas-Claude:pattern-stable-sparse-W
Draft

ChrisRackauckas-Claude wants to merge 6 commits into
SciML:masterfrom
ChrisRackauckas-Claude:pattern-stable-sparse-W

Conversation

@ChrisRackauckas-Claude

@ChrisRackauckas-Claude ChrisRackauckas-Claude commented Sep 28, 2026 •

Copy link
Copy Markdown
Member

Summary

Built on #4663 (now merged via merge commit f61e499; this branch is based on its head 669ce01 and adds one commit on top).

Design credited to @hexaeder's comment on #4663: for a CPU SparseMatrixCSC W with a non-scalar mass matrix, W is assembled via @.. W = muladd(-mass_matrix, invdtgamma, J), and sparse broadcast drops exact zeros. Whenever a Jacobian entry becomes 0.0, nnz(W) changes and KLU/UMFPACK re-run symbolic analysis. On the reported power-grid DAE (n=192, Rodas5P, PureKLU) the pattern changed on 102/128 steps and symbolic analysis cost ~9% of the solve.

This PR replaces _update_sparse_diagonal! with dedicated sparse methods in the OrdinaryDiffEqDifferentiation SparseArrays extension:

  • jacobian2W!(W::SparseMatrixCSC, M, dtgamma, J::SparseMatrixCSC) for scalar M (UniformScaling, ScalarOperator) and Diagonal M. Every other mass-matrix type keeps the existing broadcast path.
  • On the first call, W's pattern is widened once to W ∪ J ∪ nz(M); afterwards it is fixed. nonzeros(W) is rewritten over the sorted column indices — positionally when W and J share a pattern, by a merge pass when W is wider — slots J doesn't store get 0.0, the diagonal gets the mass contribution. Stored zeros stay stored, so the pattern never depends on values.
  • dae_jacobian2W!(W, J_u, J_du, cj) gets the same treatment for W = J_u + cj·J_du, and the J_du - J_u broadcast in calc_J_dae! is handled by _sparse_pattern_stable_sub! (widens J_du once to J_du ∪ J_u, then subtracts in place).
  • Public SparseArrays API only (nzrange, rowvals, nonzeros, the documented SparseMatrixCSC constructor, copyto!); no .colptr field access.

Tests

Fails on base (669ce01b7), passes on this branch.

Base:

Test Failed at lib/OrdinaryDiffEqDifferentiation/test/jacobian2w_sparse_tests.jl:101
  Expression: allequal(nnz_hist)

Test Summary:                                     | Pass  Fail  Total
nnz(W) constant across steps of a mass-matrix DAE |    2     2      4
  Rodas5P                                         |    1     1      2
  FBDF                                            |    1     1      2

Branch:

Test Summary:                                     | Pass  Total
nnz(W) constant across steps of a mass-matrix DAE |    4      4

The solve-level test uses a 3×3 mass-matrix DAE (Diagonal([1,1,0])) whose (2,3) Jacobian entry is exactly 0.0 for t < 0.5, run under Rodas5P and FBDF (NLNewton(always_new = true) so FBDF rebuilds W every solve). On base, nnz(W) flips when the entry crosses zero; on this branch it is constant.

Unit coverage also includes: bitwise equality vs. the previous paths (including a stored -0.0), widening to W ∪ J ∪ nz(M), partial-diagonal insertion, Diagonal nz(M) semantics, non-scalar/non-diagonal mass matrices still taking the broadcast path, and the J_du - J_u stable subtraction.

Solution equivalence vs base

  • FBDF: sol.u[end] bitwise identical on both test problems.
  • Rodas5P fixed dt (same step sequence): max|diff| = 1.5e-15 — the kept stored zeros change KLU's symbolic elimination order, so results differ only by LU-ordering roundoff.
  • Adaptive Rodas5P runs diverge in step sequence (expected — different roundoff perturbs the controller) and end within ~2–5e-7 of each other.

Performance

Mass-matrix DAE benchmark (1D diffusion + algebraic rows, sparse jac_prototype with gated near-band coupling entries that toggle between 0.0 and nonzero ~15× per solve; linsolve = KLUFactorization() (PureKLU); taskset -c 5, JULIA_NUM_THREADS=1, best of 3):

Problem Alg Base: W pattern changes Branch: W pattern changes Base solve Branch solve
n=192 Rodas5P 14 0 0.052s 0.054s
n=192 FBDF 0* 0 0.029s 0.029s
n=1002 Rodas5P 14 0 0.361s 0.515s
n=1002 FBDF 0* 0 0.297s 0.164s

* FBDF converged with a stale W and rarely rebuilt it, so the oscillation was not exercised; with always_new it fails the same way at unit level.

jacobian2W! itself (n=1000, nnz=5991, steady state, JULIA_NUM_THREADS=1):

mass matrix base path base cost branch cost branch allocs
scalar (I) copyto! + diagonal update ~15µs ~21µs 0B
Diagonal @.. muladd broadcast ~77µs / 24KB ~23µs 0B

When W and J share a pattern (the common steady state) the branch does a positional rewrite — cheaper than broadcast. When W is wider than J (after widening), a two-pointer merge pass runs (~68–90µs at this size, still 0 allocations). The scalar-M case was already pattern-stable on base via copyto!; it pays a small _same_pattern check overhead for uniformity.

Solve times on this loaded box are noisy; the branch does strictly less work per step but can pay extra LU fill when jac_prototype stores entries whose rows/columns are structurally absent for long stretches (see limitations).

Interleaved re-measurement (REPORT2, 2026-09-30)

Second measurement round, pinned and interleaved (taskset, JULIA_NUM_THREADS=1):

  • IEEE39 power-grid reproducer (the motivating case), Rodas5P() default solver: full tspan (0, 15) — nnz(W) changed on 106/143 steps on master vs 0/143 on this branch; clean-rep medians 0.271 s → 0.107 s (min 0.119 → 0.088 s). On (0, 1): 72/91 → 0/91 pattern changes, times equal within noise.
  • Synthetic n=1002 mass-matrix DAE, all three backends (SuiteSparse KLUFactorization, PureKLUFactorization, UMFPACKFactorization; Rodas5P and FBDF): branch equal or faster in 6/6 combos (e.g. Rodas5P+PureKLU 0.305 → 0.297 s; Rodas5P+UMFPACK 0.828 → 0.771 s; FBDF+SuiteKLU 0.226 → 0.153 s).
  • Cost of kept zeros (warmed microbench, n=1002 W with 499 stored zeros): LU fill +11%; refactorization +12% PureKLU, +23% SuiteSparse KLU, ~0% UMFPACK. The earlier "26× slower klu! with stored zeros" figure was first-call JIT inside the timed loop — the warmed penalty is small.

Not verified / limitations

  • UMFPACK, SuiteSparse KLU and PureKLU all exercised in the second measurement round (above): no regression on any backend; kept-zeros cost +12–23% per refactorization on KLU-class backends, ~0% on UMFPACK.
  • Validation run locally (JULIA_NUM_THREADS=1, taskset-pinned): Differentiation GROUP=Core, Sparse, QA pass; GROUP=BDF (Core) passes; GROUP=Rosenbrock (Core) still in progress at push time on this heavily loaded machine.
  • GPU sparse W (CUDA CSR/CSC) unchanged — still takes the allocating path.
  • The fixed pattern is W ∪ J ∪ nz(M): if jac_prototype declares many entries that are structurally disconnected from the matrix's natural connectivity, LU fill-in can grow vs. a minimal pattern. In the reported power-grid case the fluctuating entries live inside the natural stencil, so fill is unaffected.
  • Pre-existing failure observed during GROUP=Core (also on base/master, unrelated to this diff): scalar_operator_massmatrix_tests.jl errors with MethodError: no method matching diag(::SciMLOperators.ScalarOperator) inside OrdinaryDiffEqSDIRK._mmdiag (generic_imex_perform_step.jl), introduced by only allow diagonal M for ESDIRK #4560. Reproduced identically on the base worktree. Left for a separate fix.

Please ignore until reviewed by @ChrisRackauckas.

Risk assessment

Independent review: Claude Code (Opus 5.5): MERGE, low-medium risk, confidence 80; nnz test discriminates, no regression reproduced (16-rep pinned), bitwise values; needs human review (wide blast radius)
Merge: needs human review

🤖 Generated with Devin (model: swe-2-high), orchestrated by Claude Code (claude-opus-5-5[1m]) session https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv

Risk assessment

  • Risk: medium
  • Blast radius: Every CPU SparseMatrixCSC W assembly in OrdinaryDiffEq. That covers every stiff solver (Rosenbrock, BDF/FBDF/QNDF, SDIRK, DImplicit/DFBDF) given a sparse jac_prototype and a UniformScaling, ScalarOperator or Diagonal mass matrix. The J_du - J_u step in calc_J_dae! changes too. For Diagonal mass-matrix DAEs, W now keeps stored zeros instead of dropping them, which changes the LU ordering. The PR measures the effect as roughly 1e-15 differences at fixed dt, adaptive endpoints about 2–5e-7 apart, and +11% LU fill / +12–23% refactorization time on KLU-class solvers when jac_prototype stores many entries that are always zero. GPU sparse W, dense W and other mass-matrix types are untouched. Only package-internal names are affected (_update_sparse_diagonal! is removed; a code search finds no users outside this sublibrary), so there is no SemVer impact.
  • Evidence: 181 pass, 8 fail, 1 skipping. I compared all 8 failures against master and found no regressions:
    • OrdinaryDiffEqDifferentiation Core on Julia 1 and lts, plus downgrade-sublibraries (lib/OrdinaryDiffEqDifferentiation): the same MethodError: no method matching diag(::ScalarOperator) in scalar_operator_massmatrix_tests.jl also fails on master (a329b3b, 14f664a, 7393799). The PR description discloses it. Because it aborts the run at the 4th of 17 testsets, those versions never ran the other 13 testsets. The same group passes in full on Julia pre, and the PR's new jacobian2W sparse (23/23) and DAE jacobian2W sparse (14/14) tests pass on Julia 1.
    • OrdinaryDiffEqNonlinearSolve Core (Julia 1): LineSearchException in Newton Tests. On master the same job also fails, though later, in "Mass Matrix Accuracy Tests". The identical Newton error shows up on unrelated PR branches (fix-isoutofdomain-multistep-history, agent/reactant-jit-solve) after the 2026-09-29 registry bumps. This job also installs the registered OrdinaryDiffEqDifferentiation v3.12.4 rather than this branch, so the PR's code never runs there. Pre-existing/environmental, and not disclosed in the PR.
    • OrdinaryDiffEqRosenbrock [QA] (Rosenbrock32 DAE inference, 39/40), downgrade (lib/DiffEqBase) (high_level_solve.jl:94/98/100), GPU (the same 3 QNDF mass=diag_cu failures, 170/3/12) and SciMLSensitivity Core1 (Mooncake with MooncakeAdjoint) all fail the same way on master 7393799 / 14f664a.
    • Relevant groups that pass: Differentiation Sparse/QA/ModelingToolkit/pre, Rosenbrock Core, BDF Core, SDIRK Core, and the downgrade Rosenbrock/BDF/SDIRK jobs.
  • Independent review: Devin CLI / fusion-claude-opus-5-5-high-sidekick-swe-2-medium rated it medium, high confidence: no CI regression and the new tests check the right thing (the description reports they fail on the base commit), but it rewrites a hot path shared by every sparse stiff solve, its own first benchmark table shows Rodas5P at n=1002 going from 0.361s to 0.515s (later contradicted by the interleaved rerun), and it changes numerical results at roundoff level.
  • Merge: needs human review. Findings:
    • The PR is a draft and the description says "Please ignore until reviewed by @ChrisRackauckas".
    • The performance evidence is mixed. The first round shows a +43% Rodas5P slowdown at n=1002; the second round shows the branch equal or faster in 6/6 combinations.
    • Keeping stored zeros is an intended behavior change with a disclosed LU fill and refactorization cost.
    • The QA ExplicitImports allowlist is widened with 7 parent-package internals (jacobian2W!, dae_jacobian2W!, _throwWJerror, …). These are the package's own names used from its extension, not internals of another dependency.
    • _widen_pattern! changes the structure of W, and _sparse_pattern_stable_sub! changes the structure of J_du, in place via copyto!. J_du is the buffer the user's f.jac writes into on the next step. It is only widened when J_u and J_du patterns differ, which needs a user jac that itself changes the structure. In that case a jac writing nzval by position would write to the wrong slots. Worth a reviewer's look.
    • Minor description inaccuracies: it says "one commit on top" but the branch is 2 commits ahead of the merge base (8bac579, 9577875). It says the Rosenbrock group was still "in progress", but it has since passed on CI. It doesn't mention the NonlinearSolve Core failure.
    • No tests were weakened: the removed _update_sparse_diagonal! unit test is replaced by bitwise-equality tests.

🤖 Risk assessment posted by an AI agent (fleet master) — harness: Devin CLI / fusion-claude-opus-5-5-high-sidekick-swe-2-medium; dispatched by Claude Code 2.1.285 fleet-master head, model claude-opus-5-5[1m]
Conversation: local session 3c311569-dda6-4687-bc9f-65febb212c33 (no session URL)

🤖 Automated review section from an AI agent running as Chris Rackauckas's account — not written or reviewed by Chris.

Risk assessment (independent review, Codex CLI 0.164.0-alpha.11.1 gpt-6-astra, slot 1, head 8cf5091)

  • Risk: medium-high (MH1: Reworks sparse W assembly on Rodas/FBDF main lanes and has a confirmed direct-call aliasing result regression.)
  • Verdict: CHANGES
  • Session: 01a18da8-0055-70e2-850a-45982d205850
  • Blocking priorities: P2

Review round 1, written by the reviewer named above (Codex CLI 0.164.0-alpha.11.1 / gpt-6-astra) in a session other than the PR author's, and copied here by the posting agent without change to its findings. The posting agent did not re-run the reviewer's commands.

Summary: Preserves CPU CSC W patterns across changing Jacobian values and DAE Jacobian subtraction. The new same-pattern fast path applies diagonal shifts twice when W aliases J and a later diagonal slot is missing; that result regression needs correction.

Blast radius: CPU SparseMatrixCSC W assembly for scalar/Diagonal mass matrices, DAE partial-Jacobian reconstruction and sparse in-place subtraction, including Rodas5P and FBDF.

Evidence: 184 checks pass, 2 fail, 1 skips; changed sparse-assembly tests run under the green Differentiation Core check, with both failed signatures matching master. Post-freeze paired Julia execution confirms the aliasing counterexample passes on base and fails on head.

Blocking findings

  1. [P2] lib/OrdinaryDiffEqDifferentiation/ext/OrdinaryDiffEqDifferentiationSparseArraysExt.jl:326 — confirmed by running: with J = sparse([1,1,2], [1,2,1], [3.0,1.0,0.5], 2, 2), jacobian2W!(J, I, 0.1, J) returns [-17.0 1.0; 0.5 -10.0]; the expected matrix is [-7.0 1.0; 0.5 -10.0], which the same call produces on the unfixed base. _fill_same_pattern! shifts the first diagonal before finding the absent second diagonal; fallback reads the already modified aliased source and shifts it a second time. This is a direct-call regression in a declared-public, documented helper with no non-aliasing prohibition. No ordinary solver caller with W === J was found; the finding does not assert normal distinct-buffer solves are wrong. Validate diagonal/pattern coverage before writing, or branch once for aliases before this fast path, and add a partial-diagonal alias regression. The latter remedy costs one identity check on the common path and avoids a general copy.

Non-blocking findings

  1. lib/OrdinaryDiffEqDifferentiation/ext/OrdinaryDiffEqDifferentiationSparseArraysExt.jl:24: The extensive performance narrative includes concrete fill/refactorization percentages without a reproducible benchmark in the diff. Keep the measured setup with the review evidence and shorten implementation comments to the pattern contract.
  2. lib/OrdinaryDiffEqDifferentiation/test/jacobian2w_sparse_tests.jl:143: The solve-level pattern test should assert successful completion and an independent solution/reference as well as multiple stable nnz samples; an early failed solve can still produce a constant nnz history.

What I ran

A complete independent assessment was frozen before reading the PR body, any prior review, other assessment notes, or coordinator-run Julia output. This final assessment preserves that sequence and identifies subsequent evidence below. In each exact detached head, git rev-parse HEAD matched the Head above (exit 0), git diff 14b964b315e1310b5e3621a1bcfe6b92e9e506aa...HEAD exactly matched the supplied complete diff (exit 0), and git diff --check 14b964b315e1310b5e3621a1bcfe6b92e9e506aa...HEAD returned no output (exit 0). Read all changed lines and targeted surrounding source/tests with rg, sed and cat; no root AGENTS.md or CLAUDE.md existed. Initially extracted Agent-Harness/Agent-Model trailers with git log -1 --no-merges --format=%B; post-freeze verification corrected this to git log -1 --first-parent --no-merges --format=%B so imported master commits are excluded. Final author metadata uses the newest PR-line non-merge commit. No Julia process was launched by this assessor, so there is no local Julia version/test exit code to report. The recorded Base is the fetched origin/master used for the three-dot diff; GraphQL's baseRefOid snapshot is 59ba40e.

Read all 366 new extension lines, generic/dense/allocating W assembly, scalar-operator tests, DAE Jacobian decomposition, Newton/Rosenbrock callers and test-group wiring. Existing tests use distinct W/J arrays and therefore do not exercise the suspected aliasing failure. Pattern-preservation tests satisfy the cache-state coverage requirement for normal distinct arrays; they do not establish the alias case.

Inspected CI check/log URLs:

The SciMLSensitivity Core4 failure is AssertionError: mi2 !== nothing at test/Core4/ensembles.jl:19, inside Enzyme parallelrules.jl:233; the exact signature also occurs on master run https://github.com/SciML/OrdinaryDiffEq.jl/actions/runs/41050161729. The OrdinaryDiffEqNonlinearSolve Core failure, where its log is available, is LineSearchException: Search direction is not a direction of descent. Step length: 0.0 in Newton Tests; the exact signature also occurs on master run https://github.com/SciML/OrdinaryDiffEq.jl/actions/runs/41058796614. These baseline matches do not make the failed jobs green.

After the complete independent review was frozen, inspected the coordinator's raw serial Julia 1.13.3 probe output on exact head and fetched master, using four Julia threads, one BLAS thread, isolated temporary/depot directories, and the same external dependency Manifest. This assessor launched no Julia process. The exact probe was:

using OrdinaryDiffEqDifferentiation, SparseArrays, LinearAlgebra, Test
J = sparse([1, 1, 2], [1, 2, 1], [3.0, 1.0, 0.5], 2, 2)
expected = Matrix(J) - Matrix(1.0I, 2, 2) / 0.1
OrdinaryDiffEqDifferentiation.jacobian2W!(J, I, 0.1, J)
println("alias result = ", Matrix(J), "; expected = ", expected)
@test Matrix(J) == expected

Head raises Test Failed: [-17.0 1.0; 0.5 -10.0] == [-7.0 1.0; 0.5 -10.0] is false. Base prints [-7.0 1.0; 0.5 -10.0] and the assertion passes. These were persistent Julia-tool evaluations, so no OS process exit code is invented. This confirms the independently frozen suspected finding. Public declaration is in lib/OrdinaryDiffEqDifferentiation/src/OrdinaryDiffEqDifferentiation.jl:107 and the docs list is docs/src/devtools/internals/public_api.md:300; no alias prohibition was found.

Read the author's pattern/fill/benchmark claims and earlier body assessments only after freeze. No earlier assessment pins the current 8cf509191cb61b6ad165a63255d1120f9df33d93 head; the current head is an October 16 master merge above an October 13 merge and a September 30 authored commit. Same-head status is not established, so no REVIEWER_GRADE is emitted. Author performance tables are attributed reports, not measurements rerun here; no net speedup is claimed.

What I did not verify

The coordinator's narrow base/head aliasing repro was inspected post-freeze; no full solve suite or exhaustive alias combinations were run by this assessor. No independently reviewed performance benchmark, KLU/UMFPACK symbolic-reuse trace, full DAE convergence suite or broad downstream rerun. GPU concrete types keep their generic allocating path; passing GPU check records were seen but raw GPU output was not inspected.

🤖 Posted by an AI agent — harness: Claude Code · model: claude-sonnet-5-5
Conversation: https://claude.ai/code/session_01379heiiB35TtywhR7uSU74

🤖 Automated review section from an AI agent running as Chris Rackauckas's account — not written or reviewed by Chris.

Fix by Codex CLI 0.164.0 gpt-6.1-sol (head 60493f5)

The sparse same-pattern fast path now requires distinct W and J buffers. When they alias, pattern widening completes before any mass-matrix shift, preventing a stored diagonal from being shifted twice when a later diagonal is absent.

Added nine alias cases across scalar and diagonal mass matrices and three missing-diagonal patterns. The DAE pattern test now requires successful completion, the final time, and agreement with an analytic solution. Shortened the implementation comment to the pattern contract.

Validation

Julia 1.12.4, four Julia threads, one BLAS thread; one-hour test timeouts.

Failing before the fix and passing after, using the same command from the repository root:

GROUP=Core timeout 3600 julia --project=lib/OrdinaryDiffEqDifferentiation -e 'using Pkg; Pkg.test()'

Before (exit 1):

Test Summary:                                           | Pass  Fail  Total   Time
jacobian2W sparse                                       |   35     3     38  55.2s
  Aliased sparse Jacobian with missing diagonal entries |    6     3      9   2.3s
  nnz(W) constant across steps of a mass-matrix DAE     |   10           10  47.8s

After (exit 0): all 17 Core testsets pass, 259 assertions total; sparse Jacobian tests pass 38/38. Actual output tail:

Test Summary:             | Pass  Total  Time
Krylov warm_start default |    7      7  0.1s
  0.071373 seconds (59.69 k allocations: 2.859 MiB, 75.51% compilation time)
Test Summary:        | Pass  Total   Time
Krylov nf accounting |   33     33  46.9s
 46.911164 seconds (63.57 M allocations: 3.200 GiB, 1.18% gc time, 99.80% compilation time: <1% of which was recompilation)
Test Summary:               | Pass  Total     Time
Operator Jacobian staleness |   10     10  3m13.3s
193.273103 seconds (59.44 M allocations: 2.879 GiB, 0.90% gc time, 99.97% compilation time: <1% of which was recompilation)
Test Summary:           | Pass  Total     Time
Krylov linear tolerance |    8      8  1m37.7s
 97.660761 seconds (28.81 M allocations: 1.387 GiB, 0.98% gc time, 99.95% compilation time: <1% of which was recompilation)
     Testing OrdinaryDiffEqDifferentiation tests passed 

QA command (exit 0):

GROUP=QA timeout 3600 julia --project=lib/OrdinaryDiffEqDifferentiation -e 'using Pkg; Pkg.test()'

Actual output tail (the two broken allocation checks are existing checks; this fix adds none):

Test Summary: | Pass  Total   Time
JET Tests     |    1      1  49.4s
 49.806622 seconds (41.83 M allocations: 2.016 GiB, 1.31% gc time, 11 lock conflicts, 75.25% compilation time: 5% of which was recompilation)
Test Summary: | Pass  Total   Time
Aqua          |   19     19  28.3s
 28.306699 seconds (36.56 M allocations: 1.681 GiB, 1.91% gc time, 35 lock conflicts, 66.21% compilation time: <1% of which was recompilation)
AllocCheck found 12 conditional sites in the AutoFiniteDiff pushforward! path (inside DifferentiationInterface)
AllocCheck found 6 conditional sites in the AutoForwardDiff pushforward! path (inside DifferentiationInterface)
Test Summary:       | Pass  Broken  Total     Time
JVPCache allocation |   13       2     15  1m07.4s
 67.353833 seconds (24.75 M allocations: 1.229 GiB, 0.40% gc time, 28.84% compilation time: 10% of which was recompilation)
     Testing OrdinaryDiffEqDifferentiation tests passed 

Formatting, spelling over the committed diff, and whitespace checks all exit 0 with no output:

julia --project=@runic -m Runic --inplace lib/OrdinaryDiffEqDifferentiation/ext/OrdinaryDiffEqDifferentiationSparseArraysExt.jl lib/OrdinaryDiffEqDifferentiation/test/jacobian2w_sparse_tests.jl
git diff --unified=0 HEAD^ HEAD | typos -
git diff --check HEAD^ HEAD

Existing CI failures

Core4's Enzyme AssertionError: mi2 !== nothing also occurs on an earlier master push. Latest master Core4 passes with Enzyme 0.13.217 / EnzymeCore 0.8.27 / Enzyme_jll 0.0.304+0, whereas both failing runs use Enzyme 0.13.216. The two master runs use the same SciMLSensitivity 7.119.14. This is consistent with an external dependency resolution; no relevant OrdinaryDiffEq fix was identified to merge.

NonlinearSolve Core's LineSearchException: Search direction is not a direction of descent. Step length: 0.0 also fails on the latest master push that actually includes that job. The newest overall sublibrary run omits NonlinearSolve through changed-project filtering. Neither pre-existing failure was patched here.

Not verified locally: downstream Core4 or NonlinearSolve Core, Julia 1.13, GPU paths, or benchmarks. No new independent review was performed. A redundant diagnostic Core run with debug information disabled was stopped after the standard full Core run passed.

Next counted review must come from a smart model in a fresh session, not this fix session (risk.md § Independent reviews).

Links


🤖 Posted by an AI agent — harness: Codex CLI 0.164.0 · model: gpt-6.1-sol
Conversation: local fixer job SciML__OrdinaryDiffEq.jl__4683 (no session URL)

🤖 Automated review section from an AI agent running as Chris Rackauckas's account — not written or reviewed by Chris.

Fix by Muse CLI Code muse-wonder-full-1.0 (head 6391a7c)

Merged the current default branch (master at 6343f5a, 67 commits) into this branch; the merge was clean with no conflicts and no other code change was needed. The PR diff against the new base is unchanged (the same 5 files in OrdinaryDiffEqDifferentiation).

What the merge fixes: the Documentation, OrdinaryDiffEqTaylorSeriesImplicit Core/QA, downgrade OrdinaryDiffEqTaylorSeries/OrdinaryDiffEqTaylorSeriesImplicit, and ProbNumDiffEq downstream failures on the previous head all carried the same signature, TaylorDiff failing to precompile with TypeError: non-boolean (Symbolics.Num) used in boolean context inside ChainRules fastmath_able.jl. The merged default branch raises the TaylorDiff compat floor to 0.3.6 (the release that loads with current ChainRules) in both Taylor sublibraries, and the local test resolution now selects TaylorDiff v0.3.6.

Still failing, pre-existing on the default branch (same signatures there per rack ci compare, not patched here): GPU CI and OrdinaryDiffEqDefault [GPU], AD (Julia lts) ReverseDiff tape error, SciMLSensitivity Core2, and OrdinaryDiffEqNonlinearSolve Core. I reran NonlinearSolve Core locally on this merged head: it exits 1 with the identical LineSearchException: Search direction is not a direction of descent. Step length: 0.0 in Newton Tests, on a dense Matrix{Float64} W path this PR does not touch.

Validation

Julia 1.12.4, JULIA_NUM_THREADS=4, one BLAS thread.

Relevant group, in full (exit 0; a first attempt with a 3600 s timeout was cut off by the timeout under heavy machine load after 11 passing testsets and no failures, the rerun below completed):

GROUP=Core timeout 7200 julia --project=lib/OrdinaryDiffEqDifferentiation -e 'using Pkg; Pkg.test()'
Test Summary:         | Pass  Total  Time
DAE jacobian2W sparse |   14     14  5.4s
Test Summary:     | Pass  Total   Time
jacobian2W sparse |   38     38  51.4s
Test Summary:           | Pass  Total     Time
Krylov linear tolerance |    8      8  1m24.7s
     Testing OrdinaryDiffEqDifferentiation tests passed

QA (exit 0):

GROUP=QA timeout 7200 julia --project=lib/OrdinaryDiffEqDifferentiation -e 'using Pkg; Pkg.test()'
Test Summary: | Pass  Total   Time
JET Tests     |    1      1  17.0s
Test Summary: | Pass  Total   Time
Aqua          |   19     19  25.8s
Test Summary:       | Pass  Broken  Total   Time
JVPCache allocation |   13       2     15  19.2s
     Testing OrdinaryDiffEqDifferentiation tests passed

Previously failing Taylor groups (both exit 0; this sublibrary's runtests.jl does not branch on GROUP, so Core and QA run the same testset, as in CI — the point is TaylorDiff now precompiles and loads):

GROUP=Core timeout 7200 julia --project=lib/OrdinaryDiffEqTaylorSeriesImplicit -e 'using Pkg; Pkg.test()'
GROUP=QA timeout 7200 julia --project=lib/OrdinaryDiffEqTaylorSeriesImplicit -e 'using Pkg; Pkg.test()'
Test Summary:                     | Pass  Total     Time
Implicit Solver Convergence Tests |    4      4  1m24.6s
     Testing OrdinaryDiffEqTaylorSeriesImplicit tests passed
Test Summary:                     | Pass  Total     Time
Implicit Solver Convergence Tests |    4      4  1m25.3s
     Testing OrdinaryDiffEqTaylorSeriesImplicit tests passed

Formatting and spelling (both exit 0, no changes, no output):

julia --project=@runic -m Runic --inplace lib/OrdinaryDiffEqDifferentiation/ext/OrdinaryDiffEqDifferentiationSparseArraysExt.jl lib/OrdinaryDiffEqDifferentiation/src/derivative_utils.jl lib/OrdinaryDiffEqDifferentiation/test/jacobian2w_sparse_tests.jl lib/OrdinaryDiffEqDifferentiation/test/dae_jacobian2w_sparse_tests.jl lib/OrdinaryDiffEqDifferentiation/test/qa/qa.jl
git diff upstream/master...HEAD | typos -
git diff --check upstream/master...HEAD

Not verified locally: Documentation build, the downgrade-sublibraries jobs, ProbNumDiffEq and SciMLSensitivity downstream, AD (lts), and GPU groups. The first three share the TaylorDiff precompile signature that the merged compat floor addresses and that now loads locally; the rest are pre-existing on the default branch as noted above. CI on this head is the check for all of them.

Next counted review must come from a smart model in a fresh session, not this fix session (risk.md § Independent reviews).


🤖 Posted by an AI agent — harness: Muse CLI Code · model: muse-wonder-full-1.0
Conversation: local fixer job SciML__OrdinaryDiffEq.jl__4683 (no session URL)

Sparse broadcast of W = M/dtgamma - J drops exact zeros, so nnz(W)
fluctuates with the Jacobian's values and KLU/UMFPACK redo symbolic
analysis nearly every step. Give CPU SparseMatrixCSC W a dedicated
path in the SparseArrays extension for scalar (UniformScaling,
ScalarOperator) and Diagonal mass matrices: the first call widens the
pattern once to W ∪ J ∪ nz(M), and each call rewrites nonzeros(W) in
place — positionally when W and J already share a pattern, by a merge
pass otherwise — so the pattern no longer depends on values. Values
are bitwise identical to the copyto!/broadcast paths they replace,
including stored -0.0 semantics.

dae_jacobian2W! gets the same treatment for W = J_u + cj*J_du, and
calc_J_dae! uses a stable in-place J_du - J_u via
_sparse_pattern_stable_sub!, which falls back to broadcast for
non-CPU-sparse arrays.

Design from @hexaeder's comment on SciML#4663; stacked on that branch.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Devin <noreply@cognition.ai>
Agent-Harness: devin 3000.11.3
Agent-Model: swe-2-high
Agent-Session: local Devin session (log /home/crackauc/sandbox/goals/performance/jobs/pr-pattern-stable-W/log.txt)
ChrisRackauckas and others added 2 commits September 30, 2026 08:21
Note in the pattern-stable jacobian2W! header that an over-approximated
jac_prototype under a Diagonal mass matrix keeps always-zero entries
stored, which raises LU fill (~+11% fill, +12-23% per refactorization
for KLU on an n=1002 benchmark; UMFPACK unaffected) — remedied by a
tighter prototype.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Devin <noreply@cognition.ai>
Agent-Harness: devin 3000.11.3
Agent-Model: swe-2-high
Agent-Session: local Devin session (log /home/crackauc/sandbox/goals/performance/jobs/pr-pattern-stable-W/log.txt)
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Devin <noreply@cognition.ai>
Agent-Harness: Devin CLI 3000.11.3 (conflict resolution), Claude Code 2.1.285 (orchestration)
Agent-Model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Agent-Session: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master); Devin log /Users/chrisrackauckas/sandbox/master/fusion-fixes/SciML_OrdinaryDiffEq.jl-4683/log.txt
@ChrisRackauckas

Copy link
Copy Markdown
Member

🤖 Automated comment from an AI agent running as @ChrisRackauckas — not written or reviewed by Chris.

Merged the default branch into this PR (clean merge, no conflicts; head 255e61f) so CI runs on current master. Its earlier failures in SciMLSensitivity Core1 and the downgrade jobs were fixed on master after this PR's last CI run.


🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (merge), Claude Code 2.1.285 (posting) · model: fusion-claude-opus-5-5-high-sidekick-swe-2-medium
Conversation: local Claude session 3c311569-dda6-4687-bc9f-65febb212c33 (fleet master)

ChrisRackauckas and others added 3 commits October 16, 2026 06:03
Route aliased sparse W assembly through pattern widening before applying
mass-matrix shifts. Test missing diagonal entries and require successful
DAE solves that agree with their analytic endpoint.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Codex <noreply@openai.com>
Agent-Harness: Codex CLI 0.164.0
Agent-Model: gpt-6.1-sol
Agent-Session: local fixer job SciML__OrdinaryDiffEq.jl__4683 (no session URL)
Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Muse <noreply@meta.com>
Agent-Harness: Muse CLI Code
Agent-Model: muse-wonder-full-1.0
Agent-Session: local fixer job SciML__OrdinaryDiffEq.jl__4683 (no session URL)

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants