Homotopy drivers: read the interior corrector solution through typed barriers (drop last_sol union boxing)

[ChrisRackauckas-Claude]

Part of the homotopy allocation-reduction effort; gets stepping toward zero.

Root causes

Two independent per-step allocation sources were profiled for the continuation drivers (HomotopySweep, ArcLengthContinuation), measured with the SLOPE method inside a function (n = 50 coupled cubic, adaptive = false fixed-step, nsteps/maxsteps 50 vs 250, NewtonRaphson inner, per-step = (allocated@250 − allocated@50)/200):

(1) last_sol union boxing (this PR). The drivers build one inner-solver cache and re-drive it every continuation step. In the interior accept/quality block, last_sol is union-typed: it is assigned from both the concrete interior solve!(cache) and the fresh solve(inner_prob, inner) of the exempt anchor/landing solves (_sweep_exempt_solve / _arclength_fixed_solve), whose two NonlinearSolution members differ in type parameters. Reading the converged iterate as a bare last_sol.u — four times per step in the sweep (the correction/disp/scale norms and _sweep_accept!), once per step in arclength (xnew = last_sol.u) — boxed a fresh SciMLBase getproperty on every step (solution_interface.jl getproperty, ~384 B/step sweep, ~96 B/step arclength), because the two-member union is not split at the non-inlined norm_op/_sweep_accept! / chord-arithmetic call sites the value flows into.

The fix reads the needed fields (.u, .stats.nsteps) once through @noinline type-specializing field barriers (_sweep_sol_u/_sweep_sol_nsteps, _arclength_sol_u/_arclength_sol_nsteps) so the union-split happens at the barrier call (a two-member union is always split) and each specialization returns a concrete value the downstream non-inlined calls consume allocation-free. A plain same-scope solu = last_sol.u local does not suffice — it stays inferred as the union .u and still boxes downstream (this is why the read-once collapse in the closed #1042 did not remove the box; #1044 fixed the returned-solution reads at the source but not these interior-loop reads).

(2) LinearSolve refactorization ipiv (out of scope here — upstream LinearSolve, Julia 1.10 only). On Julia 1.10 each interior lu! allocates a fresh ipiv pivot vector (~900 B/step) because LinearSolve's in-place ipiv-reuse path (_reusable_lu_cacheval / _lu_reusing_ipiv!, factorization.jl) is VERSION >= v"1.11"-gated — the preallocated-ipiv LAPACK.getrf!(A, ipiv) overload only exists in Julia 1.11+. On Julia 1.11 this is already 0 B/step; standalone NewtonRaphson and a direct reinit!/solve! loop both allocate 0 there. This is a LinearSolve.jl concern (no in-tree NonlinearSolve lever) and is reported separately; it is the entire ~900 B/step gap between the Julia-1.10 and Julia-1.11 numbers below.

Before / after (per-step)

driver (NewtonRaphson inner)	Julia 1.10.11 before	after	Julia 1.11.9 before	after
HomotopySweep	1629.28 B	1245.28 B (−384.0)	96.0 B	96.0 B
ArcLengthContinuation (secant)	816.0 B	720.0 B (−96.0)	320.0 B	224.0 B (−96.0)
ArcLengthContinuation (tangent)	1344.0 B	1248.0 B (−96.0)	352.0 B	256.0 B (−96.0)

Default polyalgorithm inner also improves (the same interior union reads boxed there): sweep 4901.76 → 2917.76 B/step and arclength secant 2988.0 → 2492.0 B/step on Julia 1.10.

The −384.0 (sweep) / −96.0 (arclength) deltas are exactly the count of boxed getproperty reads removed (four vs one per step), confirmed by Profile.Allocs attribution to SciMLBase .../solution_interface.jl getproperty at the sweep lines (correction/disp/scale/accept) and the arclength xnew = last_sol.u line: 384.0 → 0.77 B/step (sweep, now one-time) and 96.0 → 0.0 B/step (arclength).

Timing (ns/step, same config): unchanged-to-improved — HomotopySweep 72.4 → 70.4 µs/step and ArcLengthContinuation secant 48.1 → 47.8 µs/step on Julia 1.10; no regression on 1.11.

Irreducible floor after this PR

On Julia 1.11 the remaining per-step allocation is 96 B/step (sweep) / 224 B/step (arclength secant), attributed by Profile.Allocs to a single site: last_sol = solve!(cache) (the interior corrector solve). The cache is union-typed because init on a fully-concrete NonlinearProblem{true, FullSpecialize, …} with NewtonRaphson has a return type that widens to a two-member {FullSpecialize, AutoSpecialize} cache union (verified: Base.return_types(init, …) returns a 2-member union even with no kwargs and a concrete FullSpecialize input). Calling solve! on that union boxes the returned NonlinearSolution (whose two members differ in size) once per step. This is a NonlinearSolveFirstOrder/SciMLBase init return-type-inference limitation upstream of these drivers — the driver already constructs a concrete FullSpecialize problem (per #1044) and does everything right; there is no driver-level lever that removes it without changing original/failure-path semantics. On Julia 1.10 this same box (~96 B) plus LinearSolve's ~900 B ipiv (cause 2) make up the residual 1245 B/step.

Behavior

Byte-identical: the adaptive sweep and arclength both still land on u = ones(50) with max error 6.66e-16 (identical to master). The failure-path reads (last_sol.resid/.retcode, original = last_sol) are left as-is — they are O(1) per solve, not per interior step, and are already gated concretely by #1044's store_original.

Tests

• All existing test/Core/homotopy_sweep_tests*.jl + arclength_tests*.jl + *_jac_tests*.jl pass unmodified: 273/273 on Julia 1.10.11 and 273/273 on Julia 1.11.9.
• Tightened arclength_tests__item8.jl: the default-inner bounds drop from 128 KB/160 KB (a pre-cache-driver artifact) to 6 KB/7 KB, and a NewtonRaphson-inner guard (< 1.1 KB secant / < 1.7 KB tangent) directly bounds the union-boxing regression.
• Added homotopy_sweep_tests__item23.jl: the sweep equivalent (default-inner < 7 KB, NewtonRaphson-inner < 1.5 KB, plus the land-on-ones correctness check).

Both regression files pass on Julia 1.10 and 1.11. Runic --check is clean on all changed files.

Note: This PR should be ignored until reviewed by @ChrisRackauckas.

🤖 Generated with Claude Code