Add warm_start option to KrylovJL GMRES/FGMRES

[ChrisRackauckas-Claude]

Important

This PR should be ignored until reviewed by @ChrisRackauckas.

Closes #70.

What

Adds a warm_start option to the Krylov.jl wrapper algorithms, currently applied to GMRES and FGMRES (other methods accept the option but ignore it):

KrylovJL_GMRES(warm_start = :hegedus)  # :none (default) | :previous | :hegedus

When the same LinearCache is solved repeatedly (cache.b = newb; solve!(cache)), the previous solution cache.u is used as the initial guess via Krylov.warm_start!:

• :previous — raw previous solution.
• :hegedus — previous solution rescaled by the Hegedüs trick x₀ = ξu, ξ = ⟨Au, b⟩/‖Au‖², which minimizes the initial residual along the previous solution's direction and guarantees ‖b - Ax₀‖ ≤ ‖b‖ (never worse than a cold start, up to the cost of one extra operator application).

Tolerance semantics are preserved: Krylov.jl measures rtol against the warm-start residual ‖M(b - Ax₀)‖ instead of ‖Mb‖, which would silently change what reltol means. The implementation folds rtol·‖Mb‖ into atol (and zeroes rtol) for warm-started solves so the stopping threshold is identical to the cold start's.

Guards: the warm start is skipped (silent cold start) for zero or nonfinite previous solutions, non-vector u (batched RHS), and non-GMRES/FGMRES workspaces.

Why / benchmarks

Newton-Krylov in OrdinaryDiffEq/NonlinearSolve currently cold-starts every GMRES solve from zero even though cache.u holds the previous increment. Benchmarks on the 2D Brusselator (full investigation with scripts and data: https://github.com/ChrisRackauckas/InternalJunk/pull/71):

• :hegedus with ILU-preconditioned Newton-Krylov: 19–44% fewer GMRES iterations across KenCarp47/FBDF, reltol 1e-3/1e-6, N=32/64. End-to-end with this branch (N=32, reltol=1e-6, ILU): KenCarp47 wall time 1.06s → 0.75s, nf 1130 → 914; FBDF nf 919 → 909.
• :previous is consistently worse than cold start for Newton increments (+8–30% iterations) — successive increments shrink in magnitude, so the raw previous solution overshoots. It is included because it is the right choice when successive solutions are close (e.g. sequences of parameterized solves, projection steps), but :hegedus is the recommended mode for Newton-type loops and :none stays the default.
• Unpreconditioned solves (~80–100 iters/solve): :hegedus is neutral (±1%).

No downstream changes are needed: OrdinaryDiffEq's dolinsolve and NonlinearSolve's linear-solve caches already carry the previous solution in cache.u, so FBDF(linsolve = KrylovJL_GMRES(warm_start = :hegedus)) works as-is (verified end-to-end against this branch with locally-widened compat).

Notes

• New public API (documented kwarg on documented algs) → minor bump to 5.1.0.
• Allocation note: :hegedus allocates one vector per solve for Au (plus ξ.*u), and preconditioned warm starts allocate one vector for ‖Mb‖. A follow-up could add scratch buffers; kept simple here for correctness. A Krylov.jl API accepting a precomputed A*x₀ would remove the duplicated matvec inside Krylov.warm_start!-based solves and halve the per-solve overhead, which matters when a strong preconditioner leaves only ~2-3 iterations per solve.
• Tests: test/Core/warm_start.jl — validation, repeated-solve correctness for all modes on GMRES/FGMRES (including operator changes), warm-start engagement (identical re-solve takes 0 iterations), the Hegedüs scaling-invariance property, preconditioned solves, complex systems, batched-RHS fallback, and zero-u fallback. Local run: 31/31 pass; GROUP=Core run locally before pushing.

Process log

• Investigated issue Better initial guesses for reinitialized iterative solvers #70; benchmarked :none/:previous/:hegedus on Brusselator via a patched LinearSolve v4.3.0 with a global toggle (see the InternalJunk PR for scripts + full tables), then implemented the winning design on main as an opt-in per-algorithm option.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

CI triage (automated)

Failing checks on this PR, all verified as pre-existing on unmodified main (not introduced by this diff):

Check	Failure	Evidence
tests / QA (1)	JET failure at test/qa/jet.jl:136 (Sparse Factorizations), 27 passed/1 failed/7 broken	identical signature+counts in main run 29449653521 at base commit 0496f13
tests / AD (lts, 1)	Static Arrays: Error During Test, 3 passed/1 errored	identical in the same main run
tests / Trim (1)	Unsatisfiable requirements ... SciMLBase resolver error	identical in the same main run
NonlinearSolve.jl/Core downstream	23 Test Problems: Broyden — problem 1/alg #4 fails norm=4.4>1e-3, problem 8/alg #4 Unexpected Pass	reproduced locally on Julia 1.10.11 with unmodified LinearSolve main + registered NonlinearSolve v4.21.1, identical values — dependency drift since main's IntegrationTest passed on 2026-07-15

Everything else is green so far, including GROUP=Core on all Julia versions (contains the new test/Core/warm_start.jl), Documentation, Downgrade, GPU, and the OrdinaryDiffEq/ModelingToolkit/SciMLSensitivity/BoundaryValueDiffEq downstream jobs.

Background investigations (bisect + root-cause) of the pre-existing main failures are running separately and will be reported when complete.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Follow-up: all 5 pre-existing failures root-caused; fixes exist

Background investigation of the failures this PR inherits from main is complete. All reproduce locally (none are runner artifacts):

Job	Bisected cause	Fix
AD (lts, 1)	#1083 (5.0 lightweight solutions): the two Failure returns of solve(::StaticLinearProblem, alg) still pass prob into the solution while Success returns pass nothing → Union return type breaks @inferred	#1094 (already open, CI green on AD; verified locally: GROUP=AD 23 Pass/1 Broken, exit 0 on 1.10.11 + 1.12.6)
QA (1)	#1083 flipped JET.@test_opt solve(prob_sparse, KLUFactorization()) from broken to plain, but the path has 24 runtime dispatches (SparseArrays checking-constructor error path + KLU getproperty)	#1095 (already open, real code fix, CI green on QA; verified locally: GROUP=QA JET 28 Pass/7 Broken, exit 0)
Trim (1)	#1071 raised the root SciMLBase floor to 3.x; test/Trim/Project.toml still capped SciMLBase at 2 → unsatisfiable	#1093 (already open, CI green on Trim; verified locally with equivalent compat: 9 Pass/1 Broken, exit 0)
NonlinearSolve.jl/Core downstream	NonlinearSolve#1039 (QN 1.14.0) replaced the triangular/LU-inverse fast paths in Broyden true-Jacobian init with a generic solve-vs-identity; machine-eps J⁻¹ difference flips two chaotic robustness tests (both versions diverge to residual ~1e179 on prob 1 — the old pass was luck). LinearSolve bit-identical 4.3.0 vs main → exonerated. The 07-15 "passing" downstream run was actually failing behind continue-on-error.	Draft PR to NonlinearSolve.jl restoring the triangular fast paths in linsolve_identity!! is in preparation (verified to restore both tests without broken-list changes)

Merging #1093, #1094, #1095 (disjoint files) should turn main's Tests workflow green, after which this PR's checks should be fully green on a rebase/rerun.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Downstream fix is now up: SciML/NonlinearSolve.jl#1091 (draft) restores the pre-NonlinearSolve#1039 exact-inverse fast paths in linsolve_identity!! for strided dense matrices. Verified locally on Julia 1.10.11: the 23-test-problems Broyden suite returns to its declared expectations with no broken-list changes (95 Pass / 20 Broken / 0 Fail), and NonlinearSolveBase Core+QA and NonlinearSolveQuasiNewton Core all pass.

Complete unblock path for this PR's CI: merge LinearSolve #1093, #1094, #1095 + NonlinearSolve #1091, then rerun checks here.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Post-rebase CI: 58/59 green

Rebased onto main with #1093/#1094/#1095 — AD, QA, and Trim now pass. The single remaining failure is NonlinearSolve.jl/Core downstream, and it is a newly unmasked pre-existing issue, again unrelated to this diff:

NonlinearSolve#1055 (2026-07-15) changed NonlinearSolveBaseForwardDiffExt (LinearSolve-based implicit sensitivities) and added the singular-Hessian NLLS Hessian #445 testset without bumping lib/NonlinearSolveBase's version (master and registry both say 2.35.0, but the registered tree lacks the ext change). The registered ext's H \ b hits Diagonal \ → SingularException(2) at forward_ad_tests__item2.jl:141-142. It was latent since 07-15 because the downstream suite previously aborted earlier at the Broyden failure (fixed by NonlinearSolve#1051); dependency versions are identical between the failing runs. LinearSolve does not appear in the failing stacktrace.

A sweep of the full NonlinearSolve master Core suite against the registered stack is running to enumerate any further latent test-vs-release mismatches, followed by a version-bump PR to NonlinearSolve.jl covering every lib whose unreleased changes master's tests depend on. Once those libs are released, this job should go green with no changes here.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Final piece: full sweep of NonlinearSolve master's Core suite (73 files) against the registered stack found exactly one latent mismatch — the known forward_ad_tests__item2.jl SingularException — and verified master's lib/NonlinearSolveBase fixes it (10/10 pass). Unreleased-changes audit: only NonlinearSolveBase (#1055 ext + #1051 descent fixes) and NonlinearSolveFirstOrder (#1051 reinit!/forcing fixes) differ from their registered trees; all other libs byte-identical.

Unblock: merge SciML/NonlinearSolve.jl#1093 (two-line version bumps: NLSBase 2.35.0→2.35.1, FirstOrder 2.2.1→2.2.2) and register; the downstream job here then resolves NLSBase 2.35.1 automatically and this PR should be 59/59.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Expanded benchmark study (4 PDEs × 4 solvers × 2 tols × 3 modes, 150 cells)

Full data, scripts, and analysis: https://github.com/ChrisRackauckas/InternalJunk/pull/71 (gmres-warmstart-issue70/benchmarks/). Problems: Brusselator 2D (N=32/64), Allen-Cahn 2D (64²), viscous Burgers 1D (1024), linear advection-diffusion 2D (64², moving source); solvers KenCarp47/TRBDF2/FBDF/Rodas5P; ILU + matrix-free JFNK configs; NonlinearSolve Newton/PseudoTransient.

Across 37 Newton-Krylov cells: :hegedus median −17% GMRES iterations (range −69%…+28%, only 3 regressions); :previous median +15% (up to +265%). Wall time: wins where solves do real Krylov work (KenCarp47+Brusselator, 5–8 iters/solve: −21…−28% time), loses where ILU is near-exact (1–3 iters/solve: median +50% time despite −14% iterations — the 2 extra matvecs + preconditioner apply per solve dominate). Unpreconditioned JFNK: neutral.

Rosenbrock hazard (new finding): Rodas5P + :previous enters a step-rejection death spiral — 0-iteration warm solves return barely-within-tolerance stage values, the error estimator rejects, dt shrinks, repeat. Observed 500× slowdown with a wrong final solution (retcode Success) on advection-diffusion. :hegedus on Rodas5P is erratic (one −39% time win, one 7.6× loss). Root cause: W-methods have no outer Newton iteration to absorb within-tolerance solve differences.

Pushed a docs commit tempering the guidance accordingly (iteration reduction reliable; time win specifically for ≳5 iters/solve Newton-Krylov; warning box against Rosenbrock use). Defaults unchanged (:none).

Follow-ups worth considering: Krylov.jl API for precomputed A*x₀ (halves overhead, would make :hegedus neutral-or-better even under strong preconditioning); SciMLBenchmarks entry; separate look at KenCarp47 loose-tol error control on advdiff/Allen-Cahn (cold-start baseline showed rel err 15 with retcode Success — unrelated to this PR).

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

API change: warm_start is now an EnumX enum, not a Symbol

Per review, replaced the Symbol-valued option with an exported, documented EnumX enum matching the existing OperatorCondition/NonstructuralZeros pattern:

KrylovJL_GMRES(warm_start = WarmStart.Hegedus)   # WarmStart.None (default) | .Previous | .Hegedus

WarmStart is exported, each variant has a docstring, and all four (WarmStart, .None, .Previous, .Hegedus) are rendered in docs/src/solvers/solvers.md. Passing a bare Symbol now errors at construction (compile-time-checked value). All three modes remain available; WarmStart.Previous's docstring flags that it is rarely appropriate for Newton loops.

Default stays WarmStart.None. I considered making Hegedus the default but the benchmark data argues against it as a global default: it is a net wall-time loss under strong preconditioning (≲3 GMRES iters/solve: median +50% time), and it is hazardous inside Rosenbrock/W-methods. A default silently changes every downstream Newton-Krylov and Rosenbrock solve, including the ones it hurts; None is the non-surprising default with Hegedus documented as the recommended opt-in. The default can be revisited once a SciMLBenchmarks entry confirms a broad win — safer to opt in than to un-break a regression.

Tests updated to the enum API (34 assertions, pass locally); Runic-clean.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

Rebased onto main (5.0.1 → 5.1.0) and marked ready for review.

Added WarmStart.Auto as the new default to support integrator-aware defaulting downstream. In standalone LinearSolve, Auto behaves exactly as None (cold start), so merging this changes no behavior on its own. A context-aware caller can resolve Auto: the companion OrdinaryDiffEq PR (in progress) resolves it to WarmStart.Hegedus for Newton-based integrators and leaves it a cold start for Rosenbrock/W-methods, without overriding a user's explicit None/Previous/Hegedus. This mirrors the existing NonstructuralZeros.Auto sentinel already in the package.

Final API: WarmStart.Auto (default) | .None | .Previous | .Hegedus, all exported and documented. Note: I don't have merge rights on this repo from this account — ready for @ChrisRackauckas to merge.

🤖 Generated with Claude Code

[ChrisRackauckas-Claude]

CI follow-up: two failing checks after the EnumX/docs commits.

• QA — genuinely caused by this PR: ExplicitImports (via SciMLTesting) can't analyze the EnumX-generated WarmStart module, exactly as it can't analyze the existing OperatorCondition/NonstructuralZeros/DefaultAlgorithmChoice enums. Fixed by adding LinearSolve.WarmStart to allow_unanalyzable in test/qa/qa.jl. Verified locally: full qa.jl now passes (Quality Assurance 18 pass / 2 pre-existing broken; Public API documentation 2/2).
• Downgrade — pre-existing on main (the umbrella Downgrade workflow is red across recent master commits d3f0790/7e3a84a0; Downgrade Sublibraries is green). This PR adds no dependencies and only a version bump, so it's unrelated.

With the QA fix pushed, this PR should be green except the pre-existing umbrella Downgrade.

🤖 Generated with Claude Code