Default Newton-Krylov integrators to Hegedus warm starting

Resolves a Krylov linear solver left at the new LinearSolve.WarmStart.Auto
default to WarmStart.Hegedus on the Newton nonlinear-solver path, so implicit
integrators using KrylovJL_GMRES()/KrylovJL_FGMRES() reuse the previous solution
(Hegedus-scaled) across the correlated linear solves of a Newton iteration.
Rosenbrock/W-method integrators never call the resolver, so their Auto solver
stays a cold start -- warm starting is unsafe there (no outer Newton iteration
absorbs the within-tolerance stage-solve perturbation, which can trigger
step-rejection feedback loops). Explicit WarmStart.None/Previous/Hegedus are
always respected.

The resolver is guarded on isdefined(LinearSolve, :WarmStart), so it is a no-op
on LinearSolve < 5.1 and this compiles unchanged there. Bumps the LinearSolve
compat of the two touched sublibraries to allow 5; umbrella-wide activation
additionally requires the OrdinaryDiffEq LinearSolve 5 compat migration across
the remaining sublibraries.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>

Author: ChrisRackauckas

lib/OrdinaryDiffEqDifferentiation/Project.toml

@@ -1,5 +1,5 @@
 name = "OrdinaryDiffEqDifferentiation"
-version = "3.4.1"
+version = "3.5.0"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
 
@@ -21,7 +21,7 @@ DiffEqDevTools = "3"
 Test = "<0.0.1, 1"
 FiniteDiff = "2.27"
 DifferentiationInterface = "0.7.18"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "3.75.0, 4, 5"
 ConstructionBase = "1.5.8"
 LinearAlgebra = "1.10"
 SciMLBase = "3.35.1"

lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl

@@ -52,6 +52,30 @@ function dolinsolve(
     return linres
 end
 
+"""
+    default_krylov_warm_start(linsolver) -> linsolver
+
+Resolve a Krylov linear solver left at `LinearSolve.WarmStart.Auto` to the mode
+appropriate for a Newton-based integrator, `LinearSolve.WarmStart.Hegedus`
+(Hegedüs-scaled reuse of the previous solution — the recommended mode for the
+sequence of correlated preconditioned solves inside a Newton iteration).
+
+Only the `Auto` default is resolved: an explicit `WarmStart.None`/`Previous`/
+`Hegedus`, a non-Krylov solver, or a LinearSolve too old to define `WarmStart`
+is returned unchanged. This is called on the Newton nonlinear-solver path only;
+Rosenbrock/W-method integrators never call it, so their `Auto` solver stays a
+cold start (warm starting is unsafe there — no outer Newton iteration absorbs
+the within-tolerance stage-solve perturbation).
+"""
+function default_krylov_warm_start(linsolver)
+    isdefined(LinearSolve, :WarmStart) || return linsolver
+    (
+        linsolver isa LinearSolve.KrylovJL &&
+            linsolver.warm_start === LinearSolve.WarmStart.Auto
+    ) || return linsolver
+    return SciMLBase.remake(linsolver; warm_start = LinearSolve.WarmStart.Hegedus)
+end
+
 """
     wrapprecs(linsolver, W, weight) -> linsolver
 

lib/OrdinaryDiffEqDifferentiation/test/runtests.jl

@@ -36,6 +36,7 @@ if TEST_GROUP ∉ ("QA", "Sparse", "ModelingToolkit")
     @time @safetestset "Autodiff Error Tests" include("autodiff_error_tests.jl")
     @time @safetestset "No Jac Tests" include("nojac_tests.jl")
     @time @safetestset "Stale W Linear Operator Tests" include("stale_w_linear_operator_tests.jl")
+    @time @safetestset "Krylov warm_start default" include("warm_start_default_tests.jl")
 end
 
 # Run sparse tests (separate environment due to ComponentArrays dep conflicts)

lib/OrdinaryDiffEqDifferentiation/test/warm_start_default_tests.jl

@@ -0,0 +1,30 @@
+using OrdinaryDiffEqDifferentiation: default_krylov_warm_start
+using LinearSolve
+using Test
+
+@testset "default_krylov_warm_start resolves Auto for the Newton path" begin
+    # Non-Krylov solvers are always returned unchanged.
+    @test default_krylov_warm_start(LUFactorization()) isa LUFactorization
+
+    if isdefined(LinearSolve, :WarmStart)
+        # Auto (the KrylovJL default) resolves to Hegedus on the Newton path.
+        @test default_krylov_warm_start(KrylovJL_GMRES()).warm_start ===
+            LinearSolve.WarmStart.Hegedus
+        # Explicit choices are respected, never overridden.
+        for ws in (
+                LinearSolve.WarmStart.None, LinearSolve.WarmStart.Previous,
+                LinearSolve.WarmStart.Hegedus,
+            )
+            @test default_krylov_warm_start(KrylovJL_GMRES(warm_start = ws)).warm_start === ws
+        end
+        # Other fields survive the remake.
+        precs = (A, p) -> (I, I)
+        resolved = default_krylov_warm_start(KrylovJL_GMRES(; precs))
+        @test resolved.warm_start === LinearSolve.WarmStart.Hegedus
+        @test resolved.precs === precs
+    else
+        # LinearSolve too old to define WarmStart: strictly a no-op.
+        k = KrylovJL_GMRES()
+        @test default_krylov_warm_start(k) === k
+    end
+end

lib/OrdinaryDiffEqNonlinearSolve/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqNonlinearSolve"
 uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>", "Yingbo Ma <mayingbo5@gmail.com>"]
-version = "2.4.0"
+version = "2.5.0"
 
 [deps]
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
@@ -58,7 +58,7 @@ FastBroadcast = "1.3"
 Random = "<0.0.1, 1"
 DiffEqDevTools = "3"
 MuladdMacro = "0.2.4"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "3.75.0, 4, 5"
 LineSearches = "7.5.1"
 LinearAlgebra = "1.10"
 OrdinaryDiffEqDifferentiation = "3"

lib/OrdinaryDiffEqNonlinearSolve/src/OrdinaryDiffEqNonlinearSolve.jl

@@ -55,7 +55,7 @@ import OrdinaryDiffEqCore: _initialize_dae!,
     @SciMLMessage
 
 import OrdinaryDiffEqDifferentiation: update_W!, is_always_new, build_uf, build_J_W,
-    WOperator, StaticWOperator, wrapprecs,
+    WOperator, StaticWOperator, wrapprecs, default_krylov_warm_start,
     build_jac_config, dolinsolve,
     resize_jac_config!, jacobian2W!, jacobian!
 

lib/OrdinaryDiffEqNonlinearSolve/src/utils.jl

@@ -432,7 +432,10 @@ function build_nlsolver(
                     )
                 end
             end
-            inner_alg = _nlalg_with_linsolve(nlalg.alg, wrapprecs(alg.linsolve, W, weight))
+            inner_alg = _nlalg_with_linsolve(
+                nlalg.alg,
+                wrapprecs(default_krylov_warm_start(alg.linsolve), W, weight)
+            )
             cache = init(prob, inner_alg, verbose = verbose.nonlinear_verbosity)
             nlcache = NonlinearSolveCache(
                 ustep, tstep, k, atmp, invγdt, prob, cache,
@@ -448,7 +451,7 @@ function build_nlsolver(
             du = isdae ? k : nothing # k will be overwritten at solve time, but has the right type.
             linprob = LinearProblem(W, _vec(k), (du, u, p, t); u0 = _vec(dz))
             linsolve = init(
-                linprob, wrapprecs(alg.linsolve, W, weight),
+                linprob, wrapprecs(default_krylov_warm_start(alg.linsolve), W, weight),
                 alias = LinearAliasSpecifier(alias_A = true, alias_b = true),
                 assumptions = LinearSolve.OperatorAssumptions(true),
                 verbose = verbose.linear_verbosity