From ef8c58fa0264d7f69f2a7e7dfb96109714bad480 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Mon, 17 Nov 2025 21:11:50 +0100
Subject: [PATCH 01/19] Eisenstat-Walker Newton-Krylov solver

---
 .../src/NonlinearSolveFirstOrder.jl           |   2 +
 .../src/eisenstat_walker.jl                   | 115 ++++++++++++++++++
 lib/NonlinearSolveFirstOrder/src/solve.jl     |  32 ++++-
 .../test/rootfind_tests.jl                    |  49 ++++++++
 src/NonlinearSolve.jl                         |   2 +-
 5 files changed, 195 insertions(+), 5 deletions(-)
 create mode 100644 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl

diff --git a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
index 72ce1ae1e..37d64f4f3 100644
--- a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
+++ b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
@@ -33,6 +33,7 @@ using ForwardDiff: ForwardDiff, Dual  # Default Forward Mode AD
 
 include("solve.jl")
 include("raphson.jl")
+include("eisenstat_walker.jl")
 include("gauss_newton.jl")
 include("levenberg_marquardt.jl")
 include("trust_region.jl")
@@ -99,6 +100,7 @@ end
 @reexport using SciMLBase, NonlinearSolveBase
 
 export NewtonRaphson, PseudoTransient
+export EisenstatWalkerNewtonKrylov
 export GaussNewton, LevenbergMarquardt, TrustRegion
 
 export RadiusUpdateSchemes
diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
new file mode 100644
index 000000000..8d767d1bd
--- /dev/null
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -0,0 +1,115 @@
+"""
+    EisenstatWalkerForcing2(; η₀, γ, α)
+
+Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described by formula (2.6).
+"""
+@concrete struct EisenstatWalkerForcing2
+    η₀
+    γ
+    α
+    safeguard
+end
+
+function EisenstatWalkerForcing2(; η₀ = 0.99, γ = 0.9, α = 2, safeguard = true)
+    EisenstatWalkerForcing2(η₀, γ, α, safeguard)
+end
+
+
+@concrete mutable struct EisenstatWalkerForcing2Cache
+    p::EisenstatWalkerForcing2
+    η
+    rnorm
+    rnorm_prev
+    internalnorm
+end
+
+
+
+function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::NonlinearSolveBase.NewtonDescentCache, J, u, fu, iter)
+    # On the first iteration we initialize η with the default initial value and stop.
+    if iter == 0
+        cache.η = cache.p.η₀
+        LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
+        return nothing
+    end
+
+    # Store previous
+    ηprev = cache.η
+
+    # Formula (2.6)
+    # ||r|| > 0 should be guaranteed by the convergence criterion
+    (; rnorm, rnorm_prev) = cache
+    (; α, γ) = cache.p
+    cache.η = γ * (rnorm / rnorm_prev)^α
+
+    # Safeguard 2 to prevent over-solving
+    if cache.p.safeguard
+        ηsg = γ*ηprev^α
+        if ηsg > 0.1 && ηsg > cache.η
+            cache.η = ηsg
+        end
+    end
+
+    # Far away from the root we also need to respect η ∈ [0,1)
+    cache.η = clamp(cache.η, 0.0, 1-eps(cache.η))
+
+    # Communicate new relative tolerance to linear solve
+    LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
+
+    return nothing
+end
+
+
+
+function post_step_forcing!(cache::EisenstatWalkerForcing2Cache, J, u, fu, δu, iter)
+    # Cache previous residual norm
+    cache.rnorm_prev = cache.rnorm
+    # And update the current one
+    cache.rnorm = cache.internalnorm(fu)
+end
+
+
+
+function InternalAPI.init(
+        prob::AbstractNonlinearProblem, alg::EisenstatWalkerForcing2, f, fu, u, p,
+        args...; internalnorm::F = L2_NORM, kwargs...
+) where {F}
+    fu_norm = internalnorm(fu)
+
+    return EisenstatWalkerForcing2Cache(
+        alg, alg.η₀, fu_norm, fu_norm, internalnorm
+    )
+end
+
+
+
+function InternalAPI.reinit!(
+        cache::EisenstatWalkerForcing2Cache; p = cache.p, kwargs...
+)
+    cache.p = p
+end
+
+
+
+"""
+    EisenstatWalkerNewtonKrylov(;
+        concrete_jac = nothing, linsolve = nothing, linesearch = missing,
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+    )
+
+An advanced Newton-Krylov implementation with support for efficient handling of sparse
+matrices via colored automatic differentiation and preconditioned linear solvers. Designed
+for large-scale and numerically-difficult nonlinear systems.
+"""
+function EisenstatWalkerNewtonKrylov(;
+        concrete_jac = nothing, linsolve::LinearSolve.AbstractKrylovSubspaceMethod, linesearch = nothing,
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+)
+    return GeneralizedFirstOrderAlgorithm(;
+        linesearch,
+        descent = NewtonDescent(; linsolve),
+        autodiff, vjp_autodiff, jvp_autodiff,
+        concrete_jac,
+        name = :EisenstatWalkerNewtonKrylov
+    )
+end
diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index 7b7267aac..df6e396a5 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -25,6 +25,7 @@ order of convergence.
     linesearch
     trustregion
     descent
+    forcing
     max_shrink_times::Int
 
     autodiff
@@ -38,12 +39,12 @@ end
 function GeneralizedFirstOrderAlgorithm(;
         descent, linesearch = missing, trustregion = missing, autodiff = nothing,
         vjp_autodiff = nothing, jvp_autodiff = nothing, max_shrink_times::Int = typemax(Int),
-        concrete_jac = Val(false), name::Symbol = :unknown
+        concrete_jac = Val(false), forcing = nothing, name::Symbol = :unknown
 )
     concrete_jac = concrete_jac isa Bool ? Val(concrete_jac) :
                    (concrete_jac isa Val ? concrete_jac : Val(concrete_jac !== nothing))
     return GeneralizedFirstOrderAlgorithm(
-        linesearch, trustregion, descent, max_shrink_times,
+        linesearch, trustregion, descent, forcing, max_shrink_times,
         autodiff, vjp_autodiff, jvp_autodiff,
         concrete_jac, name
     )
@@ -62,6 +63,7 @@ end
     # Internal Caches
     jac_cache
     descent_cache
+    forcing_cache
     linesearch_cache
     trustregion_cache
 
@@ -125,7 +127,7 @@ function InternalAPI.reinit_self!(
 end
 
 NonlinearSolveBase.@internal_caches(GeneralizedFirstOrderAlgorithmCache,
-    :jac_cache, :descent_cache, :linesearch_cache, :trustregion_cache)
+    :jac_cache, :descent_cache, :linesearch_cache, :trustregion_cache, :forcing_cache)
 
 function SciMLBase.__init(
         prob::AbstractNonlinearProblem, alg::GeneralizedFirstOrderAlgorithm, args...;
@@ -196,6 +198,7 @@ function SciMLBase.__init(
 
         has_linesearch = alg.linesearch !== missing && alg.linesearch !== nothing
         has_trustregion = alg.trustregion !== missing && alg.trustregion !== nothing
+        has_forcing = alg.forcing !== missing && alg.forcing !== nothing
 
         if has_trustregion && has_linesearch
             error("TrustRegion and LineSearch methods are algorithmically incompatible.")
@@ -204,6 +207,7 @@ function SciMLBase.__init(
         globalization = Val(:None)
         linesearch_cache = nothing
         trustregion_cache = nothing
+        forcing_cache = nothing
 
         if has_trustregion
             NonlinearSolveBase.supports_trust_region(alg.descent) ||
@@ -228,13 +232,23 @@ function SciMLBase.__init(
             globalization = Val(:LineSearch)
         end
 
+        if has_forcing
+            forcing_cache = CommonSolve.init(
+                prob, alg.forcing, fu, u; stats, internalnorm,
+                autodiff = ifelse(
+                    provided_jvp_autodiff, alg.jvp_autodiff, alg.vjp_autodiff
+                ),
+                kwargs...
+            )
+        end
+
         trace = NonlinearSolveBase.init_nonlinearsolve_trace(
             prob, alg, u, fu, J, du; kwargs...
         )
 
         cache = GeneralizedFirstOrderAlgorithmCache(
             fu, u, u_cache, prob.p, alg, prob, globalization,
-            jac_cache, descent_cache, linesearch_cache, trustregion_cache,
+            jac_cache, descent_cache, forcing_cache, linesearch_cache, trustregion_cache,
             stats, 0, maxiters, maxtime, alg.max_shrink_times, timer,
             0.0, true, termination_cache, trace, ReturnCode.Default, false, kwargs,
             initializealg, verbose
@@ -259,6 +273,12 @@ function InternalAPI.step!(
         end
     end
 
+    has_forcing = cache.forcing_cache !== nothing && cache.forcing_cache !== missing
+
+    if has_forcing
+        pre_step_forcing!(cache.forcing_cache, cache.descent_cache, J, cache.u, cache.fu, cache.nsteps)
+    end
+
     @static_timeit cache.timer "descent" begin
         if cache.trustregion_cache !== nothing &&
            hasfield(typeof(cache.trustregion_cache), :trust_region)
@@ -339,6 +359,10 @@ function InternalAPI.step!(
                    are (:LineSearch, :TrustRegion, :None)")
         end
         NonlinearSolveBase.check_and_update!(cache, cache.fu, cache.u, cache.u_cache)
+
+        if has_forcing
+            post_step_forcing!(cache.forcing_cache, J, cache.u, cache.fu, δu, cache.nsteps)
+        end
     else
         α = false
         cache.make_new_jacobian = false
diff --git a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
index 0d122c2d6..45077e262 100644
--- a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
+++ b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
@@ -4,6 +4,55 @@ include("../../../common/common_rootfind_testing.jl")
 
 end
 
+@testitem "Eisenstadt-Walker Newton-Krylov" setup=[CoreRootfindTesting] tags=[:core] begin
+    using LinearAlgebra, Random, LinearSolve
+    using BenchmarkTools: @ballocated
+    using StaticArrays: @SVector
+
+    u0s=([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+
+    @testset for (concrete_jac, linsolve) in (
+        (Val(false), KrylovJL_CG(; precs = nothing)),
+        (Val(false), KrylovJL_GMRES(; precs = nothing)),
+        (
+            Val(true),
+            KrylovJL_GMRES(;
+                    precs = (A,
+                    p = nothing) -> (
+                        Diagonal(randn!(similar(A, size(A, 1)))), LinearAlgebra.I
+                    )
+            ),
+        ),
+    )
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+            solver = EisenstatWalkerNewtonKrylov(; linsolve, concrete_jac)
+            sol = solve_oop(quadratic_f, u0; solver)
+            @test SciMLBase.successful_retcode(sol)
+            err = maximum(abs, quadratic_f(sol.u, 2.0))
+            @test err < 1e-9
+
+            cache = init(
+                NonlinearProblem{false}(quadratic_f, u0, 2.0), solver, abstol = 1e-9
+            )
+            @test (@ballocated solve!($cache)) < 200
+        end
+
+        @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+            solver = EisenstatWalkerNewtonKrylov(; linsolve, concrete_jac)
+
+            sol = solve_iip(quadratic_f!, u0; solver)
+            @test SciMLBase.successful_retcode(sol)
+            err = maximum(abs, quadratic_f(sol.u, 2.0))
+            @test err < 1e-9
+
+            cache = init(
+                NonlinearProblem{true}(quadratic_f!, u0, 2.0), solver, abstol = 1e-9
+            )
+            @test (@ballocated solve!($cache)) ≤ 64
+        end
+    end
+end
+
 @testitem "NewtonRaphson" setup=[CoreRootfindTesting] tags=[:core] begin
     using ADTypes, LineSearch, LinearAlgebra, Random, LinearSolve
     using LineSearches: LineSearches
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index dc6a20507..e511f83d3 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -21,7 +21,7 @@ using StaticArraysCore: StaticArray
 
 # Default Algorithm
 using NonlinearSolveFirstOrder: NewtonRaphson, TrustRegion, LevenbergMarquardt, GaussNewton,
-                                RUS, RobustMultiNewton
+                                RUS, RobustMultiNewton, EisenstadtWalkerNewtonKrylov
 using NonlinearSolveQuasiNewton: Broyden, Klement
 using SimpleNonlinearSolve: SimpleBroyden, SimpleKlement
 

From 8cbe2bec3501add787de963143dfe1899b41f57a Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 12:43:37 +0100
Subject: [PATCH 02/19] Missing default

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 8d767d1bd..284a350e4 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -103,13 +103,14 @@ for large-scale and numerically-difficult nonlinear systems.
 """
 function EisenstatWalkerNewtonKrylov(;
         concrete_jac = nothing, linsolve::LinearSolve.AbstractKrylovSubspaceMethod, linesearch = nothing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing, forcing = EisenstatWalkerForcing2(),
 )
     return GeneralizedFirstOrderAlgorithm(;
         linesearch,
         descent = NewtonDescent(; linsolve),
         autodiff, vjp_autodiff, jvp_autodiff,
         concrete_jac,
+        forcing,
         name = :EisenstatWalkerNewtonKrylov
     )
 end

From 8d0f0c4f7fb696d1b5d299f0c272df64420cc7f9 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 13:19:34 +0100
Subject: [PATCH 03/19] Fix typos

---
 lib/NonlinearSolveFirstOrder/src/solve.jl | 4 ++--
 src/NonlinearSolve.jl                     | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index df6e396a5..633103a0d 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -233,8 +233,8 @@ function SciMLBase.__init(
         end
 
         if has_forcing
-            forcing_cache = CommonSolve.init(
-                prob, alg.forcing, fu, u; stats, internalnorm,
+            forcing_cache = InternalAPI.init(
+                prob, alg.forcing, fu, u, u, prob.p; stats, internalnorm,
                 autodiff = ifelse(
                     provided_jvp_autodiff, alg.jvp_autodiff, alg.vjp_autodiff
                 ),
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index e511f83d3..b4bacef7d 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -21,7 +21,7 @@ using StaticArraysCore: StaticArray
 
 # Default Algorithm
 using NonlinearSolveFirstOrder: NewtonRaphson, TrustRegion, LevenbergMarquardt, GaussNewton,
-                                RUS, RobustMultiNewton, EisenstadtWalkerNewtonKrylov
+                                RUS, RobustMultiNewton, EisenstatWalkerNewtonKrylov
 using NonlinearSolveQuasiNewton: Broyden, Klement
 using SimpleNonlinearSolve: SimpleBroyden, SimpleKlement
 

From c5e8aee343a1f0e6195830795a0a7d2de8def43b Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 13:27:20 +0100
Subject: [PATCH 04/19] Add missing dispatch

---
 .../ext/NonlinearSolveBaseLinearSolveExt.jl                   | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl b/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl
index 50512300b..ee6db6760 100644
--- a/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl
+++ b/lib/NonlinearSolveBase/ext/NonlinearSolveBaseLinearSolveExt.jl
@@ -77,4 +77,8 @@ function set_lincache_A!(lincache, new_A)
     return
 end
 
+function LinearSolve.update_tolerances!(cache::LinearSolveJLCache; kwargs...)
+    LinearSolve.update_tolerances!(cache.lincache; kwargs...)
+end
+
 end

From e315dce1d72791e925034297415808707c7d093e Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:10:31 +0100
Subject: [PATCH 05/19] Add verbosity

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 5 ++++-
 lib/NonlinearSolveFirstOrder/src/solve.jl            | 1 +
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 284a350e4..11eccb124 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -21,6 +21,7 @@ end
     rnorm
     rnorm_prev
     internalnorm
+    verbosity
 end
 
 
@@ -53,6 +54,8 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Far away from the root we also need to respect η ∈ [0,1)
     cache.η = clamp(cache.η, 0.0, 1-eps(cache.η))
 
+    @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbose, :linear_verbosity)
+
     # Communicate new relative tolerance to linear solve
     LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
 
@@ -72,7 +75,7 @@ end
 
 function InternalAPI.init(
         prob::AbstractNonlinearProblem, alg::EisenstatWalkerForcing2, f, fu, u, p,
-        args...; internalnorm::F = L2_NORM, kwargs...
+        args...; verbose, internalnorm::F = L2_NORM, kwargs...
 ) where {F}
     fu_norm = internalnorm(fu)
 
diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index 633103a0d..6344561a3 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -238,6 +238,7 @@ function SciMLBase.__init(
                 autodiff = ifelse(
                     provided_jvp_autodiff, alg.jvp_autodiff, alg.vjp_autodiff
                 ),
+                verbose,
                 kwargs...
             )
         end

From 6e0861ae15655aff3aeeaf42ccdd60f20a9fdf4c Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:15:40 +0100
Subject: [PATCH 06/19] Add maximum eta to params.

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 11eccb124..734e7fffb 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -5,12 +5,13 @@ Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described
 """
 @concrete struct EisenstatWalkerForcing2
     η₀
+    ηₘₐₓ
     γ
     α
     safeguard
 end
 
-function EisenstatWalkerForcing2(; η₀ = 0.99, γ = 0.9, α = 2, safeguard = true)
+function EisenstatWalkerForcing2(; η₀ = 0.99, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true)
     EisenstatWalkerForcing2(η₀, γ, α, safeguard)
 end
 
@@ -52,7 +53,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     end
 
     # Far away from the root we also need to respect η ∈ [0,1)
-    cache.η = clamp(cache.η, 0.0, 1-eps(cache.η))
+    cache.η = clamp(cache.η, 0.0, cache.ηₘₐₓ)
 
     @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbose, :linear_verbosity)
 

From 251f8ccf60f41644873928d907faeff42f3941f8 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:24:42 +0100
Subject: [PATCH 07/19] More parameters

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 9 +++++----
 1 file changed, 5 insertions(+), 4 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 734e7fffb..ba8b1e41d 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -9,10 +9,11 @@ Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described
     γ
     α
     safeguard
+    safeguard_threshold
 end
 
-function EisenstatWalkerForcing2(; η₀ = 0.99, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true)
-    EisenstatWalkerForcing2(η₀, γ, α, safeguard)
+function EisenstatWalkerForcing2(; η₀ = 0.99, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
+    EisenstatWalkerForcing2(η₀, ηₘₐₓ, γ, α, safeguard, safeguard_threshold)
 end
 
 
@@ -47,7 +48,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Safeguard 2 to prevent over-solving
     if cache.p.safeguard
         ηsg = γ*ηprev^α
-        if ηsg > 0.1 && ηsg > cache.η
+        if ηsg > cache.p.safeguard_threshold && ηsg > cache.η
             cache.η = ηsg
         end
     end
@@ -55,7 +56,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Far away from the root we also need to respect η ∈ [0,1)
     cache.η = clamp(cache.η, 0.0, cache.ηₘₐₓ)
 
-    @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbose, :linear_verbosity)
+    @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
 
     # Communicate new relative tolerance to linear solve
     LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)

From adc159e4cd2590c68467335a6b0019e0f5be481b Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:25:16 +0100
Subject: [PATCH 08/19] Lower initial eta

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index ba8b1e41d..88199a3a2 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -12,7 +12,7 @@ Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described
     safeguard_threshold
 end
 
-function EisenstatWalkerForcing2(; η₀ = 0.99, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
+function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
     EisenstatWalkerForcing2(η₀, ηₘₐₓ, γ, α, safeguard, safeguard_threshold)
 end
 

From 9770292a1bdc7786eae6e887c82f38a3e3994ebb Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:31:40 +0100
Subject: [PATCH 09/19] Derp

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 88199a3a2..387e8a623 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -82,7 +82,7 @@ function InternalAPI.init(
     fu_norm = internalnorm(fu)
 
     return EisenstatWalkerForcing2Cache(
-        alg, alg.η₀, fu_norm, fu_norm, internalnorm
+        alg, alg.η₀, fu_norm, fu_norm, internalnorm, verbose
     )
 end
 

From 910d1cdca8a815a2bf1b9235044d8f0720cb9449 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 14:39:05 +0100
Subject: [PATCH 10/19] Derp again

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 387e8a623..8afc5a160 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -54,7 +54,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     end
 
     # Far away from the root we also need to respect η ∈ [0,1)
-    cache.η = clamp(cache.η, 0.0, cache.ηₘₐₓ)
+    cache.η = clamp(cache.η, 0.0, cache.p.ηₘₐₓ)
 
     @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
 

From a2f3071ece6fd0f05c33bab5ffc4b695833183dc Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 16:52:01 +0100
Subject: [PATCH 11/19] Debug instability

---
 lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl    | 8 +++++---
 lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl | 3 ++-
 lib/NonlinearSolveFirstOrder/src/raphson.jl             | 7 +++++--
 lib/NonlinearSolveFirstOrder/src/solve.jl               | 8 ++++----
 lib/NonlinearSolveFirstOrder/src/trust_region.jl        | 4 ++--
 5 files changed, 18 insertions(+), 12 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 8afc5a160..7e828bc1a 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -12,7 +12,7 @@ Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described
     safeguard_threshold
 end
 
-function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.99, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
+function EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
     EisenstatWalkerForcing2(η₀, ηₘₐₓ, γ, α, safeguard, safeguard_threshold)
 end
 
@@ -32,6 +32,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # On the first iteration we initialize η with the default initial value and stop.
     if iter == 0
         cache.η = cache.p.η₀
+        @SciMLMessage("Eisenstat-Walker initial iteration to η=$(cache.η).", cache.verbosity, :linear_verbosity)
         LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
         return nothing
     end
@@ -56,7 +57,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Far away from the root we also need to respect η ∈ [0,1)
     cache.η = clamp(cache.η, 0.0, cache.p.ηₘₐₓ)
 
-    @SciMLMessage("Eisenstat-Walker update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
+    @SciMLMessage("Eisenstat-Walker iter $iter update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
 
     # Communicate new relative tolerance to linear solve
     LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
@@ -69,8 +70,9 @@ end
 function post_step_forcing!(cache::EisenstatWalkerForcing2Cache, J, u, fu, δu, iter)
     # Cache previous residual norm
     cache.rnorm_prev = cache.rnorm
-    # And update the current one
     cache.rnorm = cache.internalnorm(fu)
+
+    # @SciMLMessage("Eisenstat-Walker sanity check: $(cache.internalnorm(fu + J*δu)) ≤ $(cache.η * cache.internalnorm(fu)).", cache.verbosity, :linear_verbosity)
 end
 
 
diff --git a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
index 66906fe9e..102d50e30 100644
--- a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
+++ b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
@@ -57,7 +57,8 @@ function LevenbergMarquardt(;
         vjp_autodiff,
         jvp_autodiff,
         name = :LevenbergMarquardt,
-        concrete_jac = Val(true)
+        concrete_jac = Val(true),
+
     )
 end
 
diff --git a/lib/NonlinearSolveFirstOrder/src/raphson.jl b/lib/NonlinearSolveFirstOrder/src/raphson.jl
index d482f4846..07cfbb068 100644
--- a/lib/NonlinearSolveFirstOrder/src/raphson.jl
+++ b/lib/NonlinearSolveFirstOrder/src/raphson.jl
@@ -1,7 +1,8 @@
 """
     NewtonRaphson(;
         concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
+        forcing = nothing,
     )
 
 An advanced NewtonRaphson implementation with support for efficient handling of sparse
@@ -10,13 +11,15 @@ for large-scale and numerically-difficult nonlinear systems.
 """
 function NewtonRaphson(;
         concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
+        forcing = nothing,
 )
     return GeneralizedFirstOrderAlgorithm(;
         linesearch,
         descent = NewtonDescent(; linsolve),
         autodiff, vjp_autodiff, jvp_autodiff,
         concrete_jac,
+        forcing,
         name = :NewtonRaphson
     )
 end
diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index 6344561a3..ae8ff95c8 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -314,6 +314,10 @@ function InternalAPI.step!(
     δu, descent_intermediates = descent_result.δu, descent_result.extras
 
     if descent_result.success
+        if has_forcing
+            post_step_forcing!(cache.forcing_cache, J, cache.u, cache.fu, δu, cache.nsteps)
+        end
+
         cache.make_new_jacobian = true
         if cache.globalization isa Val{:LineSearch}
             @static_timeit cache.timer "linesearch" begin
@@ -360,10 +364,6 @@ function InternalAPI.step!(
                    are (:LineSearch, :TrustRegion, :None)")
         end
         NonlinearSolveBase.check_and_update!(cache, cache.fu, cache.u, cache.u_cache)
-
-        if has_forcing
-            post_step_forcing!(cache.forcing_cache, J, cache.u, cache.fu, δu, cache.nsteps)
-        end
     else
         α = false
         cache.make_new_jacobian = false
diff --git a/lib/NonlinearSolveFirstOrder/src/trust_region.jl b/lib/NonlinearSolveFirstOrder/src/trust_region.jl
index 03a1b3f9a..553748331 100644
--- a/lib/NonlinearSolveFirstOrder/src/trust_region.jl
+++ b/lib/NonlinearSolveFirstOrder/src/trust_region.jl
@@ -6,7 +6,7 @@
         shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
         shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
         max_shrink_times::Int = 32,
-        vjp_autodiff = nothing, autodiff = nothing, jvp_autodiff = nothing
+        vjp_autodiff = nothing, autodiff = nothing, jvp_autodiff = nothing,
     )
 
 An advanced TrustRegion implementation with support for efficient handling of sparse
@@ -29,7 +29,7 @@ function TrustRegion(;
         shrink_threshold::Real = 1 // 4, expand_threshold::Real = 3 // 4,
         shrink_factor::Real = 1 // 4, expand_factor::Real = 2 // 1,
         max_shrink_times::Int = 32,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
+        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing,
 )
     descent = Dogleg(; linsolve)
     trustregion = GenericTrustRegionScheme(;

From 3ddabe5785e8c509c75190b76cde4d6db0f8901c Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 19:59:26 +0100
Subject: [PATCH 12/19] Verbosity

---
 lib/NonlinearSolveBase/src/verbosity.jl        | 18 ++++++++++++------
 .../src/eisenstat_walker.jl                    |  6 ++++--
 2 files changed, 16 insertions(+), 8 deletions(-)

diff --git a/lib/NonlinearSolveBase/src/verbosity.jl b/lib/NonlinearSolveBase/src/verbosity.jl
index ab4b558f1..dec5f140a 100644
--- a/lib/NonlinearSolveBase/src/verbosity.jl
+++ b/lib/NonlinearSolveBase/src/verbosity.jl
@@ -68,6 +68,7 @@ verbose = NonlinearVerbosity(
     termination_condition
     # Numerical
     threshold_state
+    forcing
 end
 
 # Group classifications
@@ -76,7 +77,7 @@ const error_control_options = (
     :termination_condition
 )
 const performance_options = ()
-const numerical_options = (:threshold_state,)
+const numerical_options = (:threshold_state,:forcing)
 
 function option_group(option::Symbol)
     if option in error_control_options
@@ -138,7 +139,8 @@ function NonlinearVerbosity(;
         alias_u0_immutable = WarnLevel(),
         linsolve_failed_noncurrent = WarnLevel(),
         termination_condition = WarnLevel(),
-        threshold_state = WarnLevel()
+        threshold_state = WarnLevel(),
+        forcing = Silent(),
     )
 
     # Apply group-level settings
@@ -173,7 +175,8 @@ function NonlinearVerbosity(verbose::AbstractVerbosityPreset)
             alias_u0_immutable = Silent(),
             linsolve_failed_noncurrent = WarnLevel(),
             termination_condition = Silent(),
-            threshold_state = Silent()
+            threshold_state = Silent(),
+            forcing = Silent(),
         )
     elseif verbose isa Standard
         # Standard: Everything from Minimal + non-fatal warnings
@@ -186,7 +189,8 @@ function NonlinearVerbosity(verbose::AbstractVerbosityPreset)
             alias_u0_immutable = WarnLevel(),
             linsolve_failed_noncurrent = WarnLevel(),
             termination_condition = WarnLevel(),
-            threshold_state = WarnLevel()
+            threshold_state = WarnLevel(),
+            forcing = InfoLevel(),
         )
     elseif verbose isa All
         # All: Maximum verbosity - every possible logging message at InfoLevel
@@ -196,7 +200,8 @@ function NonlinearVerbosity(verbose::AbstractVerbosityPreset)
             alias_u0_immutable = WarnLevel(),
             linsolve_failed_noncurrent = WarnLevel(),
             termination_condition = WarnLevel(),
-            threshold_state = InfoLevel()
+            threshold_state = InfoLevel(),
+            forcing = InfoLevel(),
         )
     end
 end
@@ -208,7 +213,8 @@ end
         Silent(),
         Silent(),
         Silent(),
-        Silent()
+        Silent(),
+        Silent(),
     )
 end
 
diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 7e828bc1a..d01ec3095 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -29,10 +29,12 @@ end
 
 
 function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::NonlinearSolveBase.NewtonDescentCache, J, u, fu, iter)
+    @SciMLMessage("Eisenstat-Walker forcing residual norm $(cache.rnorm) with rate estimate $(cache.rnorm / cache.rnorm_prev).", cache.verbosity, :forcing)
+
     # On the first iteration we initialize η with the default initial value and stop.
     if iter == 0
         cache.η = cache.p.η₀
-        @SciMLMessage("Eisenstat-Walker initial iteration to η=$(cache.η).", cache.verbosity, :linear_verbosity)
+        @SciMLMessage("Eisenstat-Walker initial iteration to η=$(cache.η).", cache.verbosity, :forcing)
         LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)
         return nothing
     end
@@ -57,7 +59,7 @@ function pre_step_forcing!(cache::EisenstatWalkerForcing2Cache, descend_cache::N
     # Far away from the root we also need to respect η ∈ [0,1)
     cache.η = clamp(cache.η, 0.0, cache.p.ηₘₐₓ)
 
-    @SciMLMessage("Eisenstat-Walker iter $iter update to η=$(cache.η).", cache.verbosity, :linear_verbosity)
+    @SciMLMessage("Eisenstat-Walker iter $iter update to η=$(cache.η).", cache.verbosity, :forcing)
 
     # Communicate new relative tolerance to linear solve
     LinearSolve.update_tolerances!(descend_cache.lincache; reltol=cache.η)

From 5d026f763d8006cb234c302732ca326044aaab61 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 20:32:53 +0100
Subject: [PATCH 13/19] Export forcing

---
 lib/NonlinearSolveBase/src/NonlinearSolveBase.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
index 43d7e8439..e1922ab5d 100644
--- a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
+++ b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
@@ -95,6 +95,8 @@ export RelTerminationMode, AbsTerminationMode,
 export DescentResult, SteepestDescent, NewtonDescent, DampedNewtonDescent, Dogleg,
        GeodesicAcceleration
 
+export EisenstatWalkerForcing2
+
 export NonlinearSolvePolyAlgorithm
 
 export NonlinearVerbosity

From 5aef765d2795e563f4fb8ddeaf3d1d13ee7fec1b Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 20:38:03 +0100
Subject: [PATCH 14/19] Bump LinearSolve to reflect that the new API is used

---
 Project.toml | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/Project.toml b/Project.toml
index 95ca7785b..92b356633 100644
--- a/Project.toml
+++ b/Project.toml
@@ -24,6 +24,7 @@ PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -85,7 +86,7 @@ LeastSquaresOptim = "0.8.5"
 LineSearch = "0.1.4"
 LineSearches = "7.3"
 LinearAlgebra = "1.10"
-LinearSolve = "3.46"
+LinearSolve = "3.48"
 MINPACK = "1.2"
 MPI = "0.20.22"
 NLSolvers = "0.5"
@@ -106,9 +107,9 @@ Random = "1.10"
 ReTestItems = "1.24"
 Reexport = "1.2.2"
 ReverseDiff = "1.15"
-SciMLLogging = "1.3"
 SIAMFANLEquations = "1.0.1"
 SciMLBase = "2.127"
+SciMLLogging = "1.3"
 SimpleNonlinearSolve = "2.11"
 SparseArrays = "1.10"
 SparseConnectivityTracer = "1"
@@ -146,9 +147,9 @@ PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"

From 1ff888fe3426245708c428be50ef965d4684ea45 Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 20:52:38 +0100
Subject: [PATCH 15/19] Remove export of helper constructor

---
 .../src/NonlinearSolveFirstOrder.jl           |  1 -
 .../src/eisenstat_walker.jl                   | 26 -------------------
 .../test/rootfind_tests.jl                    |  4 +--
 3 files changed, 2 insertions(+), 29 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
index 37d64f4f3..4cf18354f 100644
--- a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
+++ b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
@@ -100,7 +100,6 @@ end
 @reexport using SciMLBase, NonlinearSolveBase
 
 export NewtonRaphson, PseudoTransient
-export EisenstatWalkerNewtonKrylov
 export GaussNewton, LevenbergMarquardt, TrustRegion
 
 export RadiusUpdateSchemes
diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index d01ec3095..7bf780db9 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -97,29 +97,3 @@ function InternalAPI.reinit!(
 )
     cache.p = p
 end
-
-
-
-"""
-    EisenstatWalkerNewtonKrylov(;
-        concrete_jac = nothing, linsolve = nothing, linesearch = missing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing
-    )
-
-An advanced Newton-Krylov implementation with support for efficient handling of sparse
-matrices via colored automatic differentiation and preconditioned linear solvers. Designed
-for large-scale and numerically-difficult nonlinear systems.
-"""
-function EisenstatWalkerNewtonKrylov(;
-        concrete_jac = nothing, linsolve::LinearSolve.AbstractKrylovSubspaceMethod, linesearch = nothing,
-        autodiff = nothing, vjp_autodiff = nothing, jvp_autodiff = nothing, forcing = EisenstatWalkerForcing2(),
-)
-    return GeneralizedFirstOrderAlgorithm(;
-        linesearch,
-        descent = NewtonDescent(; linsolve),
-        autodiff, vjp_autodiff, jvp_autodiff,
-        concrete_jac,
-        forcing,
-        name = :EisenstatWalkerNewtonKrylov
-    )
-end
diff --git a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
index 45077e262..f6e6d1dcf 100644
--- a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
+++ b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
@@ -25,7 +25,7 @@ end
         ),
     )
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            solver = EisenstatWalkerNewtonKrylov(; linsolve, concrete_jac)
+            solver = NewtonKrylov(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
             sol = solve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             err = maximum(abs, quadratic_f(sol.u, 2.0))
@@ -38,7 +38,7 @@ end
         end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-            solver = EisenstatWalkerNewtonKrylov(; linsolve, concrete_jac)
+            solver = NewtonKrylov(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
 
             sol = solve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)

From 06d1d1a7b04a1b99ebb441a683a45d805f6ca44e Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Tue, 18 Nov 2025 22:21:18 +0100
Subject: [PATCH 16/19] Docstring for forcing alg

---
 .../src/eisenstat_walker.jl                      | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
index 7bf780db9..66030044c 100644
--- a/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
+++ b/lib/NonlinearSolveFirstOrder/src/eisenstat_walker.jl
@@ -1,7 +1,19 @@
 """
-    EisenstatWalkerForcing2(; η₀, γ, α)
+    EisenstatWalkerForcing2(; η₀ = 0.5, ηₘₐₓ = 0.9, γ = 0.9, α = 2, safeguard = true, safeguard_threshold = 0.1)
 
-Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described by formula (2.6).
+Algorithm 2 from the classical work by Eisenstat and Walker (1996) as described by formula (2.6):
+    ηₖ = γ * (||rₖ|| / ||rₖ₋₁||)^α
+
+Here the variables denote:
+    rₖ residual at iteration k
+    η₀   ∈ [0,1) initial value for η
+    ηₘₐₓ ∈ [0,1) maximum value for η
+    γ    ∈ [0,1) correction factor
+    α    ∈ [1,2) correction exponent
+
+Furthermore, the proposed safeguard is implemented:
+    ηₖ = max(ηₖ, γ*ηₖ₋₁^α) if γ*ηₖ₋₁^α > safeguard_threshold
+to prevent ηₖ from shrinking too fast.
 """
 @concrete struct EisenstatWalkerForcing2
     η₀

From 4f06060d0f0c36f8b6d702b443bbc1dad19fb79c Mon Sep 17 00:00:00 2001
From: termi-official <termi-official@users.noreply.github.com>
Date: Wed, 26 Nov 2025 21:17:47 +0100
Subject: [PATCH 17/19] Fix CI

---
 lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl | 3 +--
 lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl     | 4 ++--
 src/NonlinearSolve.jl                                   | 2 +-
 3 files changed, 4 insertions(+), 5 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
index 102d50e30..66906fe9e 100644
--- a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
+++ b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
@@ -57,8 +57,7 @@ function LevenbergMarquardt(;
         vjp_autodiff,
         jvp_autodiff,
         name = :LevenbergMarquardt,
-        concrete_jac = Val(true),
-
+        concrete_jac = Val(true)
     )
 end
 
diff --git a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
index f6e6d1dcf..df6d6a8ec 100644
--- a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
+++ b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
@@ -25,7 +25,7 @@ end
         ),
     )
         @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
-            solver = NewtonKrylov(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
+            solver = NewtonRaphson(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
             sol = solve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)
             err = maximum(abs, quadratic_f(sol.u, 2.0))
@@ -38,7 +38,7 @@ end
         end
 
         @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
-            solver = NewtonKrylov(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
+            solver = NewtonRaphson(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
 
             sol = solve_iip(quadratic_f!, u0; solver)
             @test SciMLBase.successful_retcode(sol)
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index b4bacef7d..dc6a20507 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -21,7 +21,7 @@ using StaticArraysCore: StaticArray
 
 # Default Algorithm
 using NonlinearSolveFirstOrder: NewtonRaphson, TrustRegion, LevenbergMarquardt, GaussNewton,
-                                RUS, RobustMultiNewton, EisenstatWalkerNewtonKrylov
+                                RUS, RobustMultiNewton
 using NonlinearSolveQuasiNewton: Broyden, Klement
 using SimpleNonlinearSolve: SimpleBroyden, SimpleKlement
 

From 56b8808e5ac3eac53903ac83db5cbbf622795757 Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Thu, 27 Nov 2025 09:57:42 +0100
Subject: [PATCH 18/19] Fix some corner cases

---
 lib/NonlinearSolveFirstOrder/src/solve.jl           | 4 ++--
 lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl | 4 +---
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index ae8ff95c8..ed99c305f 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -198,7 +198,7 @@ function SciMLBase.__init(
 
         has_linesearch = alg.linesearch !== missing && alg.linesearch !== nothing
         has_trustregion = alg.trustregion !== missing && alg.trustregion !== nothing
-        has_forcing = alg.forcing !== missing && alg.forcing !== nothing
+        has_forcing = alg.forcing !== missing && alg.forcing !== nothing && !(u isa Number) && !(J isa Diagonal)
 
         if has_trustregion && has_linesearch
             error("TrustRegion and LineSearch methods are algorithmically incompatible.")
@@ -274,7 +274,7 @@ function InternalAPI.step!(
         end
     end
 
-    has_forcing = cache.forcing_cache !== nothing && cache.forcing_cache !== missing
+    has_forcing = cache.forcing_cache !== nothing && cache.forcing_cache !== missing && !(cache.u isa Number) && !(J isa Diagonal)
 
     if has_forcing
         pre_step_forcing!(cache.forcing_cache, cache.descent_cache, J, cache.u, cache.fu, cache.nsteps)
diff --git a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
index df6d6a8ec..6d4b72e87 100644
--- a/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
+++ b/lib/NonlinearSolveFirstOrder/test/rootfind_tests.jl
@@ -9,8 +9,6 @@ end
     using BenchmarkTools: @ballocated
     using StaticArrays: @SVector
 
-    u0s=([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
-
     @testset for (concrete_jac, linsolve) in (
         (Val(false), KrylovJL_CG(; precs = nothing)),
         (Val(false), KrylovJL_GMRES(; precs = nothing)),
@@ -24,7 +22,7 @@ end
             ),
         ),
     )
-        @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0])
             solver = NewtonRaphson(; forcing=EisenstatWalkerForcing2(), linsolve, concrete_jac)
             sol = solve_oop(quadratic_f, u0; solver)
             @test SciMLBase.successful_retcode(sol)

From 53b0bff3bf133bbee54585e64ed4f8c8c3eb919c Mon Sep 17 00:00:00 2001
From: termi-official <9196588+termi-official@users.noreply.github.com>
Date: Thu, 27 Nov 2025 10:03:48 +0100
Subject: [PATCH 19/19] Move export to the right place

---
 lib/NonlinearSolveBase/src/NonlinearSolveBase.jl             | 2 --
 lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl | 2 ++
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
index e1922ab5d..43d7e8439 100644
--- a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
+++ b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl
@@ -95,8 +95,6 @@ export RelTerminationMode, AbsTerminationMode,
 export DescentResult, SteepestDescent, NewtonDescent, DampedNewtonDescent, Dogleg,
        GeodesicAcceleration
 
-export EisenstatWalkerForcing2
-
 export NonlinearSolvePolyAlgorithm
 
 export NonlinearVerbosity
diff --git a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
index 4cf18354f..c2dc141db 100644
--- a/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
+++ b/lib/NonlinearSolveFirstOrder/src/NonlinearSolveFirstOrder.jl
@@ -102,6 +102,8 @@ end
 export NewtonRaphson, PseudoTransient
 export GaussNewton, LevenbergMarquardt, TrustRegion
 
+export EisenstatWalkerForcing2
+
 export RadiusUpdateSchemes
 
 export GeneralizedFirstOrderAlgorithm