[ChrisRackauckas-Claude]

Please ignore until reviewed by @ChrisRackauckas.

What changed and why

construct_jacobian_cache handed the solver's own residual array to JacobianCache as the output buffer, and (cache::JacobianCache)(u) passes that buffer as y to DI.jacobian!. For an in-place problem differentiated with ForwardDiff, every Jacobian evaluation therefore overwrote the solver's fu with the dual-number primal of f(u).

That primal is not bitwise equal to a plain f(u) call. For example, sin of a Dual goes through sincos. On a 100-dim dense test problem (du = A*u + 150 .* sin.(u) .- p), 74 of 100 entries of fu change in the last bits after jac_cache(u0). The Jacobian itself is bitwise identical. So after each Jacobian evaluation, the residual used by the Newton step and by the termination check silently became the AD primal.

This PR passes the cache's existing private fu_cache instead. fu_cache is already allocated by construct_jacobian_cache for DI preparation, so no buffer is added. The diff is one line in lib/NonlinearSolveBase/src/jacobian.jl, plus a test.

This is also the prerequisite for the follow-up PR that reuses the init-time Jacobian on the first Newton step. That PR can only be bitwise identical to its base if a Jacobian evaluation no longer has this side effect.

Behavior change: in-place AD iterates change at roundoff level

In-place problems with an AD Jacobian now step from the plain residual instead of the AD primal, so their iterates change at roundoff level. I compared master and this branch on 98 solver/problem configurations:

• Problems: dense 100 IIP and OOP, dense IIP with a user jac, sparse 1000 with a jac_prototype and with a sparsity detector, SVector{4}, scalar.
• Algorithms: NR, NR with FiniteDiff, NR + GMRES (lazy and concrete), TR, TR(Hei), LM, PseudoTransient, Broyden/Klement with true-Jacobian init, and two polyalgorithms with alias_u0 = true.
• Each configuration covers a full solve, the iterates of a step! loop, and solve! after reinit! with a new u0/p, a new p, and the original values.

87 of 98 configurations are bitwise identical in every iterate, residual, retcode and stat. That includes all OOP, SVector, scalar, sparse, user-jac and matrix-free cases. The 11 that differ are all dense in-place ForwardDiff cases:

case	retcode master/PR	nsteps	final-u rel. diff	max iterate rel. diff
NR	Success/Success	5/5	1.3e-17	4.6e-15
NR + GMRES (concrete J)	Success/Success	5/5	4.7e-17	4.7e-15
TR	Success/Success	4/4	6.8e-17	1.0e-13
TR (Hei)	Success/Success	13/13	2.1e-16	4.7e-14
LM	Success/Success	11/11	1.6e-14	2.7e-13
PseudoTransient	Success/Success	11/11	3.4e-19	1.9e-15
Broyden (true J init)	Success/Success	8/8	0	8.4e-15
Klement (true J diagonal)	Success/Success	10/10	0	3.4e-15
PolyAlg(PT, NR)	Success/Success	13/13	8.8e-17	0
PolyAlg(PT, Broyden true J)	Success/Success	16/16	8.1e-17	0
Klement (true J), diverging	MaxIters/Unstable	1000/239	—	—

The last row is a solve that diverges on both master and this branch: on master its iterates reach 1e259 before MaxIters. A roundoff-level change alters how it fails. No converging case changed retcode or step count.

Cost

The buffer already existed, so there is no new allocation. Pinned with taskset, dense 100 IIP, solve(prob, alg):

alg	init bytes master → PR	solve bytes	min µs master → PR
NewtonRaphson	183240 → 183208	263296 → 263296	3799 → 3805
TrustRegion	274248 → 274248	274248 → 274248	5710 → 5705
LevenbergMarquardt	444608 → 444608	1296904 → 1296904	12466 → 12452

Verification

New test, lib/NonlinearSolveBase/test/jacobian_private_fu.jl: after jac_cache(u), the residual array equals a plain f!(fu, u, p) call bitwise. It also checks the Jacobian against the analytic one.

Failing before the fix (origin/master 0055e314):

A: Test Failed at .../lib/NonlinearSolveBase/test/jacobian_private_fu.jl:20
  Expression: fu == fu_plain
   Evaluated: [239.26381734800142, 239.6068173480014, ...] == [239.26381734800134, 239.6068173480013, ...]
Test Summary: | Pass  Fail  Total  Time
A             |    2     2      4  7.6s

Passing after:

Test Summary: | Pass  Total  Time
A             |    4      4  4.9s

Full suites on this branch (Julia 1.12):

GROUP=Core  (lib/NonlinearSolveBase)  ... Jacobian evaluation leaves the residual untouched |    4      4  2.4s
     Testing NonlinearSolveBase tests passed
GROUP=QA    (lib/NonlinearSolveBase)  QA |   20     20  45.1s
     Testing NonlinearSolveBase tests passed

Runic (Runic --check) and typos are clean on the changed files.

Not verified

• Downstream sublibrary test suites (FirstOrder, QuasiNewton, SpectralMethods, top-level NonlinearSolve) were not run on this PR alone. Because of machine load, the FirstOrder and QuasiNewton Core groups were run on the stacked follow-up PR, which contains this change. Any test that compares in-place AD iterates bitwise against stored values would see the roundoff change.
• Julia LTS and pre, GPU, and docs build were not run.
• Only ForwardDiff was measured. Other DI backends either never write the primal into y or now write it into a private buffer, so they can only become more consistent.
• No version bump. Following the repo convention, NonlinearSolveBase is released in a separate PR.

Risk assessment

• Risk: low
• Blast radius: JacobianCache for in-place problems with an AD Jacobian (dense and sparse) in every solver built on construct_jacobian_cache. No public API change. Iterates of in-place AD solves change at roundoff level, as tabulated above.
• Evidence: the new test fails before and passes after; the 98-configuration master comparison above; NonlinearSolveBase Core and QA pass; no allocation or time change.
• Independent review: pending
• Merge: needs human review. It intentionally changes numerical iterates, even if only at roundoff.

🤖 Generated with Claude Code (model: claude-opus-5-5) https://claude.ai/code/session_01LPHREnnonfLg1VcE1EJovv

Independent review: Devin Fusion (fusion-claude-opus-5-5-high-sidekick-swe-2-medium) rated it low, verdict MERGE: #1365 (comment)

[ChrisRackauckas]

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

Independent review (Devin CLI 3000.11.3, model fusion-claude-opus-5-5-high-sidekick-swe-2-medium): MERGE, risk low.

Full review

VERDICT: MERGE
RISK: low

PR: #1365 (head 7dd4590, base origin/master 8a86577). No linked issue; stacked follow-up is #1367.

The change is one line: construct_jacobian_cache hands JacobianCache its private fu_cache instead of the caller's residual fu, so in-place DI.jacobian! stops overwriting the solver residual with the AD primal. fu_cache was already allocated and already used for DI preparation, so execution and preparation now see the same buffer. I checked every jac_cache(u) call site (FirstOrder step!, native bounded init/step, QuasiNewton TrueJacobianInitialization, arclength tangent, construct_extension_jac). None of them relies on the Jacobian call refreshing fu: each one evaluates f at the same u first (FirstOrder settles deferred residuals in refresh_residual! before the Jacobian). Line search and trust-region caches evaluate trial points into their own buffers. Running the probes below confirmed this.

Blocking findings

None.

Non-blocking findings

1. The regression test depends on a ForwardDiff roundoff quirk (lib/NonlinearSolveBase/test/jacobian_private_fu.jl:16-21). It only detects the bug because sin(::Dual)'s primal differs from sin in the last bits. Confirmed by running: it fails on master (2 of 4 fail, fu == fu_plain) and passes on the PR. But if ForwardDiff/DI ever make the primal bitwise equal, the test passes on the buggy code too. A sentinel works for any backend: fill!(fu, 12345.0); jac_cache(u); @test all(==(12345.0), fu), plus @test jac_cache.fu !== fu. My probe scratch/probe_sentinel.jl does this. It fails on master and passes on the PR.
2. The sparse path is affected but not tested. AGENTS.md asks for sparse prototype coverage. The change also covers the AutoSparse path (jac_prototype + SparseMatrixColorings). Confirmed by running: the sentinel probe passes on the PR for sparse jac_prototype with ForwardDiff and FiniteDiff, and J isa SparseMatrixCSC holds. It would cost little to loop the new test over a dense and a sparse case.
3. The blast-radius section leaves out the extension solvers. construct_extension_jac (lib/NonlinearSolveBase/src/wrappers.jl:217) passes the extension's own resid buffer, which is LSO's y and FastLM's LMWorkspace residual. Confirmed by running scratch/probe_ext.jl: on master, calling the Jacobian at [3,4] overwrote resid with f([3,4]) = [8,15]. On the PR, resid is left alone. This is the same fix and it looks harmless. NLsolve (trust/newton), LeastSquaresOptim (lm/dogleg) and FastLevenbergMarquardt give identical retcodes, step counts and residuals on master and PR; FastLM differs only at roundoff. The PR body should still list it. (MINPACK errors with cfunction: closures are not supported on this platform on both master and the PR on this aarch64 Mac. That failure is pre-existing and not caused by this PR, but MINPACK went untested here.)
4. The 98-configuration comparison table has no script or output attached. I could not rerun it as written. My independent 20-configuration subset matches its pattern (see What I ran): only dense in-place ForwardDiff cases differ, the differences are at roundoff, and retcodes and step counts are identical. My TR max per-iterate relative difference was 3.8e-11 against the body's 1.0e-13. That number is componentwise-relative and driven by near-zero entries; the final-u relative difference is 1.6e-16.
5. The "Full suites on this branch" heading overstates the evidence: only the NonlinearSolveBase Core and QA groups were run. The "Not verified" section does say this.
6. NonlinearSolveBase has no version bump. That matches the repo convention: git log -- lib/NonlinearSolveBase/Project.toml shows separate "Release NonlinearSolveBase vX" PRs.

What I ran

All runs used Julia 1.12.7, with two scratch envs built from git archive of origin/master and of the PR head. Both resolved to an identical Manifest; only the dev paths differ.

• The author's test jacobian_private_fu.jl. On master: fails, fu == fu_plain false, the same output as the PR body. On the PR: 4/4 pass.
• scratch/probe_sentinel.jl (fill fu with a sentinel, call the Jacobian, check fu is untouched). Master: fails on dense ForwardDiff. PR: passes dense FD, dense FiniteDiff, sparse-proto FD and sparse-proto FiniteDiff, and the sparse J stays SparseMatrixCSC.
• scratch/probe_solvers.jl + compare.jl, run on master and PR. Problems: dense 100-dim IIP (A*u + 150 sin u - p) and IIP Rosenbrock. Algorithms: NR, NR+BackTracking, TR, LM, PseudoTransient, Broyden(true J), Klement(true J), NR+GMRES concrete J, default polyalg. Each configuration ran a step! loop that checked cache.fu == f(cache.u) bitwise after every step: 0 mismatches on the PR. Results: the dense IIP cases differ at roundoff (final-u rel ≤ 6.4e-14, which is LM; the others are ~1.6e-16), with identical retcode and nsteps. Rosenbrock and the default polyalg are bitwise identical. (Broyden+BackTracking hit a config error from my setup on both branches, so it is excluded.)
• scratch/probe_ext.jl and scratch/probe_extsolve.jl: results as in non-blocking finding 3.
• On the PR tree: GROUP=Core Pkg.test() for NonlinearSolveQuasiNewton passed, run from a scratch copy so nothing was written into the checkout. Test files run directly also passed: lib/NonlinearSolveFirstOrder/test/native_bounded_tests.jl (6006/6006), test/Core/arclength_jac_tests__item1.jl (43/43), lib/NonlinearSolveBase/test/extension_jac.jl (10/10), lib/NonlinearSolveQuasiNewton/test/core_tests__item12.jl (reinit with true-J init, 63/63).
• Runic --check on the 3 changed files is clean. typos on the changed files is clean.

What I did not verify

• The full NonlinearSolveFirstOrder Core group, SpectralMethods, the top-level NonlinearSolve test groups, and NonlinearSolveBase QA. I ran only the targeted files above, because of host limits.
• Julia LTS/pre, GPU arrays, Enzyme/Mooncake/ReverseDiff backends, the PETSc/SIAMFANL/NLSolvers extensions, and MINPACK (which errors on this platform independently of the PR).
• The PR body's allocation and timing table.

---

🤖 Posted by an AI agent — harness: Devin CLI 3000.11.3 (review), 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)