Newton-Krylov integrators default to Hegedüs warm start (require LinearSolve 5.1) (#3991)

* 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>

* Make downgrade compat floors consistent in the two touched sublibraries

Two stale minimum-version declarations were mutually unsatisfiable, which is
what the Downgrade Sublibraries CI hits:

- LinearSolve floor 3.75.0 was inconsistent with the NonlinearSolveBase 2.34.1
  floor, which itself requires LinearSolve >= 4.3.0. Raised to 4.3.0.
- OrdinaryDiffEqNonlinearSolve declared SciMLOperators >= 1.24.3 while its path
  dependency OrdinaryDiffEqDifferentiation requires >= 1.24.4. Raised to 1.24.4.

Neither drops any actually-installable configuration. The warm-start resolver
remains guarded on isdefined(LinearSolve, :WarmStart), so it still no-ops on
LinearSolve 4.x and activates on 5.1+.

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

* Require LinearSolve 5.1 across the monorepo and drop the WarmStart guard

Makes the Newton-Krylov Hegedus warm-start default unconditionally live instead
of silently dormant:

- OrdinaryDiffEqNonlinearSolve and OrdinaryDiffEqDifferentiation require
  LinearSolve = "5.1" (where WarmStart exists), and default_krylov_warm_start
  drops its isdefined(LinearSolve, :WarmStart) guard.
- Every other monorepo sublibrary that declared a LinearSolve compat is moved to
  "5.1" as well. These all depend (transitively) on NonlinearSolve/Differentiation,
  so they can never resolve with LinearSolve < 5.1 anyway; keeping a 3/4 lower
  bound would be dead metadata and would re-break the Downgrade CI (which floors
  deps to their minimum). Test sub-environments updated to match.
- Minor version bumps for each sublibrary whose compat changed.

Verified: fresh resolution pins LinearSolve 5.1.0; resolver unit test passes
(7/7); FBDF and Rodas5P solve with a default KrylovJL_GMRES(); the resolver maps
Auto -> Hegedus with the guard removed.

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

---------

Co-authored-by: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Co-authored-by: Claude Fable 5 <noreply@anthropic.com>

Author: 3 people

lib/DelayDiffEq/Project.toml

@@ -1,7 +1,7 @@
 name = "DelayDiffEq"
 uuid = "bcd4f6db-9728-5f36-b5f7-82caef46ccdb"
 authors = ["Chris Rackauckas <accounts@chrisrackauckas.com>"]
-version = "6.0.7"
+version = "6.1.0"
 
 [deps]
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
@@ -70,7 +70,7 @@ FastBroadcast = "1.3"
 FiniteDiff = "2.27"
 ForwardDiff = "1.3.3"
 LinearAlgebra = "1"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "5.1"
 Logging = "1"
 OrdinaryDiffEqBDF = "2"
 OrdinaryDiffEqCore = "4"

lib/OrdinaryDiffEqAMF/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqAMF"
 uuid = "08082164-f6a1-4363-b3df-3aa6fcf571ad"
 authors = ["Utkarsh <rajpututkarsh530@gmail.com>"]
-version = "2.1.1"
+version = "2.2.0"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -15,7 +15,7 @@ SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SciMLTesting = "2.1"
 LinearAlgebra = "1.10"
 SafeTestsets = "0.1.0"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "5.1"
 OrdinaryDiffEqCore = "4"
 Reexport = "1.2.2"
 SciMLBase = "3.35.1"

lib/OrdinaryDiffEqBDF/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqBDF"
 uuid = "6ad6398a-0878-4a85-9266-38940aa047c8"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "2.3.2"
+version = "2.4.0"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
@@ -56,7 +56,7 @@ Random = "<0.0.1, 1"
 DiffEqDevTools = "3"
 DifferentiationInterface = "0.7.18"
 MuladdMacro = "0.2.4"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "5.1"
 PrecompileTools = "1.2.1, 1.3"
 LinearAlgebra = "1.10"
 OrdinaryDiffEqDifferentiation = "3"

lib/OrdinaryDiffEqDefault/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqDefault"
 uuid = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "2.3.0"
+version = "2.4.0"
 
 [deps]
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
@@ -50,7 +50,7 @@ Random = "<0.0.1, 1"
 DiffEqDevTools = "3"
 OrdinaryDiffEqBDF = "2"
 OrdinaryDiffEqVerner = "2"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "5.1"
 PrecompileTools = "1.2.1, 1.3"
 EnumX = "1.0.4"
 LinearAlgebra = "1.10"

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 = "5.1"
 ConstructionBase = "1.5.8"
 LinearAlgebra = "1.10"
 SciMLBase = "3.35.1"

lib/OrdinaryDiffEqDifferentiation/src/linsolve_utils.jl

@@ -52,6 +52,29 @@ 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` or a non-Krylov solver 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)
+    (
+        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

@@ -37,6 +37,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/sparse/Project.toml

@@ -13,5 +13,5 @@ OrdinaryDiffEqSDIRK = {path = "../../../OrdinaryDiffEqSDIRK"}
 
 [compat]
 ComponentArrays = "0.15, 1"
-LinearSolve = "3, 4"
+LinearSolve = "5.1"
 OrdinaryDiffEq = "7"

lib/OrdinaryDiffEqDifferentiation/test/warm_start_default_tests.jl

@@ -0,0 +1,24 @@
+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
+
+    # 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
+end

lib/OrdinaryDiffEqExponentialRK/Project.toml

@@ -1,7 +1,7 @@
 name = "OrdinaryDiffEqExponentialRK"
 uuid = "e0540318-69ee-4070-8777-9e2de6de23de"
 authors = ["ParamThakkar123 <paramthakkar864@gmail.com>"]
-version = "2.1.1"
+version = "2.2.0"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
@@ -48,7 +48,7 @@ Random = "<0.0.1, 1"
 DiffEqDevTools = "3"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqVerner = "2"
-LinearSolve = "3.75.0, 4"
+LinearSolve = "5.1"
 ExponentialUtilities = "1.27"
 LinearAlgebra = "1.10"
 OrdinaryDiffEqDifferentiation = "3"