From 3840930e65b4e0440176c52921c475da97a9e71c Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Mon, 21 Aug 2023 15:52:05 -0400
Subject: [PATCH 1/6] Add support for JFNK

---
 Project.toml                    |  3 ++-
 docs/Project.toml               |  2 ++
 docs/src/tutorials/nonlinear.md | 19 ++++++++++++++
 src/jacobian.jl                 | 45 +++++++++++++++++++++------------
 src/levenberg.jl                | 17 ++++++-------
 src/raphson.jl                  | 21 ++++++---------
 src/trustRegion.jl              | 17 ++++++-------
 src/utils.jl                    | 18 +++++--------
 test/basictests.jl              | 15 ++++++-----
 9 files changed, 90 insertions(+), 67 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8e204f043..f180547dc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -41,6 +41,7 @@ julia = "1.6"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -48,4 +49,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics"]
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve"]
diff --git a/docs/Project.toml b/docs/Project.toml
index c5c7721db..f6889de90 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,6 +1,7 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveMINPACK = "c100e077-885d-495a-a2ea-599e143bf69d"
 SciMLNLSolve = "e9a6253c-8580-4d32-9898-8661bb511710"
@@ -12,6 +13,7 @@ Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 [compat]
 BenchmarkTools = "1"
 Documenter = "0.27"
+LinearSolve = "2"
 NonlinearSolve = "1"
 NonlinearSolveMINPACK = "0.1"
 SciMLNLSolve = "0.1"
diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md
index 0d59cfc39..bf4ba3da4 100644
--- a/docs/src/tutorials/nonlinear.md
+++ b/docs/src/tutorials/nonlinear.md
@@ -31,3 +31,22 @@ uspan = (1.0, 2.0) # brackets
 probB = IntervalNonlinearProblem(f, uspan)
 sol = solve(probB, Falsi())
 ```
+
+## Using Jacobian Free Newton Krylov (JNFK) Methods
+
+If we want to solve the first example, without constructing the entire Jacobian
+
+```@example
+using NonlinearSolve, LinearSolve
+
+function f!(res, u, p)
+    @. res = u * u - p
+end
+u0 = [1.0, 1.0]
+p = 2.0
+probN = NonlinearProblem(f!, u0, p)
+
+linsolve = LinearSolve.KrylovJL_GMRES()
+sol = solve(probN, NewtonRaphson(), reltol = 1e-9, maxiters=4)
+sol = solve(probN, NewtonRaphson(; linsolve), reltol = 1e-9, maxiters=4)
+```
diff --git a/src/jacobian.jl b/src/jacobian.jl
index cb696400e..85c6479ad 100644
--- a/src/jacobian.jl
+++ b/src/jacobian.jl
@@ -24,6 +24,8 @@ function jacobian_finitediff!(J, f, x, jac_config, cache)
     2 * maximum(jac_config.colorvec))
 end
 
+# NoOp for Jacobian if it is not a Abstract Array -- For eg, JacVec Operator
+jacobian!(J, cache) = J
 function jacobian!(J::AbstractMatrix{<:Number}, cache)
     f = cache.f
     uf = cache.uf
@@ -52,14 +54,16 @@ function jacobian!(J::AbstractMatrix{<:Number}, cache)
     nothing
 end
 
-function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2}
+function build_jac_and_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2}
     haslinsolve = hasfield(typeof(alg), :linsolve)
 
-    if !SciMLBase.has_jac(f) && # No Jacobian if has analytical solution
-       ((concrete_jac(alg) === nothing && (!haslinsolve || (haslinsolve && # No Jacobian if linsolve doesn't want it
-           (alg.linsolve === nothing || LinearSolve.needs_concrete_A(alg.linsolve))))) ||
-        (concrete_jac(alg) !== nothing && concrete_jac(alg))) # Jacobian if explicitly asked for
-        jac_prototype = f.jac_prototype
+    has_analytic_jac = SciMLBase.has_jac(f)
+    linsolve_needs_jac = (concrete_jac(alg) === nothing &&
+                          (!haslinsolve || (haslinsolve && (alg.linsolve === nothing ||
+                             LinearSolve.needs_concrete_A(alg.linsolve)))))
+    alg_wants_jac = (concrete_jac(alg) !== nothing && concrete_jac(alg))
+
+    if !has_analytic_jac && (linsolve_needs_jac || alg_wants_jac)
         sparsity, colorvec = sparsity_colorvec(f, u)
 
         if alg_autodiff(alg)
@@ -70,25 +74,34 @@ function build_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F1, F2}
             else
                 typeof(ForwardDiff.Tag(uf, eltype(u)))
             end
-            jac_config = ForwardColorJacCache(uf, u, _chunksize; colorvec = colorvec,
-                sparsity = sparsity, tag = T)
+            jac_config = ForwardColorJacCache(uf, u, _chunksize; colorvec, sparsity,
+                tag = T)
         else
             if alg_difftype(alg) !== Val{:complex}
-                jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg),
-                    colorvec = colorvec,
-                    sparsity = sparsity)
+                jac_config = FiniteDiff.JacobianCache(tmp, du1, du2, alg_difftype(alg);
+                    colorvec, sparsity)
             else
                 jac_config = FiniteDiff.JacobianCache(Complex{eltype(tmp)}.(tmp),
-                    Complex{eltype(du1)}.(du1), nothing,
-                    alg_difftype(alg), eltype(u),
-                    colorvec = colorvec,
-                    sparsity = sparsity)
+                    Complex{eltype(du1)}.(du1), nothing, alg_difftype(alg), eltype(u);
+                    colorvec, sparsity)
             end
         end
     else
         jac_config = nothing
     end
-    jac_config
+
+    J = if !linsolve_needs_jac
+        # We don't need to construct the Jacobian
+        JacVec(uf, u; autodiff = AutoFiniteDiff())
+    else
+        if f.jac_prototype === nothing
+            ArrayInterface.undefmatrix(u)
+        else
+            f.jac_prototype
+        end
+    end
+
+    return J, jac_config
 end
 
 function get_chunksize(jac_config::ForwardDiff.JacobianConfig{
diff --git a/src/levenberg.jl b/src/levenberg.jl
index 92dbe5f29..f3a48f90f 100644
--- a/src/levenberg.jl
+++ b/src/levenberg.jl
@@ -228,21 +228,20 @@ end
 
 function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{true})
     uf = JacobianWrapper(f, p)
-    J = ArrayInterface.undefmatrix(u)
+
+    du1 = zero(u)
+    du2 = zero(u)
+    tmp = zero(u)
+    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
 
     linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
     weight = similar(u)
-    recursivefill!(weight, false)
+    # Q: Setting this to false leads to residual = 0 in GMRES?
+    recursivefill!(weight, true)
 
     Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
             nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    jac_config = build_jac_config(alg, f, uf, du1, u, tmp, du2)
+    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
 
     uf, linsolve, J, du1, jac_config
 end
diff --git a/src/raphson.jl b/src/raphson.jl
index 48f224959..572c68499 100644
--- a/src/raphson.jl
+++ b/src/raphson.jl
@@ -106,25 +106,20 @@ end
 
 function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true})
     uf = JacobianWrapper(f, p)
-    J = if f.jac_prototype === nothing
-        ArrayInterface.undefmatrix(u)
-    else
-        f.jac_prototype
-    end
+
+    du1 = zero(u)
+    du2 = zero(u)
+    tmp = zero(u)
+    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
 
     linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
     weight = similar(u)
-    recursivefill!(weight, false)
+    # Q: Setting this to false leads to residual = 0 in GMRES?
+    recursivefill!(weight, true)
 
     Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
             nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    jac_config = build_jac_config(alg, f, uf, du1, u, tmp, du2)
+    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
 
     uf, linsolve, J, du1, jac_config
 end
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
index d070bdfc8..7c331cf97 100644
--- a/src/trustRegion.jl
+++ b/src/trustRegion.jl
@@ -280,21 +280,20 @@ end
 
 function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true})
     uf = JacobianWrapper(f, p)
-    J = ArrayInterface.undefmatrix(u)
+
+    du1 = zero(u)
+    du2 = zero(u)
+    tmp = zero(u)
+    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
 
     linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
     weight = similar(u)
-    recursivefill!(weight, false)
+    # Q: Setting this to false leads to residual = 0 in GMRES?
+    recursivefill!(weight, true)
 
     Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
             nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true,
-        Pl = Pl, Pr = Pr)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    jac_config = build_jac_config(alg, f, uf, du1, u, tmp, du2)
+    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
 
     uf, linsolve, J, du1, jac_config
 end
diff --git a/src/utils.jl b/src/utils.jl
index 9eb7cc6f1..9d72e230f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -51,8 +51,7 @@ function alg_difftype(alg::AbstractNewtonAlgorithm{
     FDT,
     ST,
     CJ,
-}) where {CS, AD, FDT,
-    ST, CJ}
+}) where {CS, AD, FDT, ST, CJ}
     FDT
 end
 
@@ -62,8 +61,7 @@ function concrete_jac(alg::AbstractNewtonAlgorithm{
     FDT,
     ST,
     CJ,
-}) where {CS, AD, FDT,
-    ST, CJ}
+}) where {CS, AD, FDT, ST, CJ}
     CJ
 end
 
@@ -73,9 +71,7 @@ function get_chunksize(alg::AbstractNewtonAlgorithm{
     FDT,
     ST,
     CJ,
-}) where {CS, AD,
-    FDT,
-    ST, CJ}
+}) where {CS, AD, FDT, ST, CJ}
     Val(CS)
 end
 
@@ -85,17 +81,15 @@ function standardtag(alg::AbstractNewtonAlgorithm{
     FDT,
     ST,
     CJ,
-}) where {CS, AD, FDT,
-    ST, CJ}
+}) where {CS, AD, FDT, ST, CJ}
     ST
 end
 
 DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing
 
 function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing,
-    du = nothing, u = nothing, p = nothing, t = nothing,
-    weight = nothing, cachedata = nothing,
-    reltol = nothing) where {P}
+    du = nothing, u = nothing, p = nothing, t = nothing, weight = nothing,
+    cachedata = nothing, reltol = nothing) where {P}
     A !== nothing && (linsolve.A = A)
     b !== nothing && (linsolve.b = b)
     linu !== nothing && (linsolve.u = linu)
diff --git a/test/basictests.jl b/test/basictests.jl
index ea792dd4f..241faaf0f 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -46,9 +46,9 @@ sol = benchmark_scalar(sf, csu0)
 # @test (@ballocated benchmark_mutable(ff, cu0)) < 200
 # @test (@ballocated benchmark_scalar(sf, csu0)) < 400
 
-function benchmark_inplace(f, u0)
+function benchmark_inplace(f, u0, linsolve)
     probN = NonlinearProblem{true}(f, u0)
-    solver = init(probN, NewtonRaphson(), abstol = 1e-9)
+    solver = init(probN, NewtonRaphson(; linsolve), abstol = 1e-9)
     sol = solve!(solver)
 end
 
@@ -57,9 +57,11 @@ function ffiip(du, u, p)
 end
 u0 = [1.0, 1.0]
 
-sol = benchmark_inplace(ffiip, u0)
-@test sol.retcode === ReturnCode.Success
-@test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+for linsolve in (nothing, KrylovJL_GMRES())
+    sol = benchmark_inplace(ffiip, u0, linsolve)
+    @test sol.retcode === ReturnCode.Success
+    @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+end
 
 u0 = [1.0, 1.0]
 probN = NonlinearProblem{true}(ffiip, u0)
@@ -209,8 +211,7 @@ end
 
 function benchmark_inplace(f, u0, radius_update_scheme)
     probN = NonlinearProblem{true}(f, u0)
-    solver = init(probN, TrustRegion(radius_update_scheme = radius_update_scheme),
-        abstol = 1e-9)
+    solver = init(probN, TrustRegion(; radius_update_scheme), abstol = 1e-9)
     sol = solve!(solver)
 end
 

From 695840f9f222992342364b2faabe41d2f8623de2 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Mon, 21 Aug 2023 16:12:57 -0400
Subject: [PATCH 2/6] Missing LinearSolve import

---
 docs/src/solvers/BracketingSolvers.md | 2 +-
 docs/src/tutorials/nonlinear.md       | 3 +--
 test/basictests.jl                    | 1 +
 3 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/solvers/BracketingSolvers.md b/docs/src/solvers/BracketingSolvers.md
index 0aeae323b..8205fff0f 100644
--- a/docs/src/solvers/BracketingSolvers.md
+++ b/docs/src/solvers/BracketingSolvers.md
@@ -7,7 +7,7 @@ Solves for ``f(t)=0`` in the problem defined by `prob` using the algorithm
 
 ## Recommended Methods
 
-`ITP()` is the recommended method for the scalar interval root-finding problems. It is particularly well-suited for cases where the function is smooth and well-behaved; and achieved superlinear convergence while retaining the optimal worst-case performance of the Bisection method. For more details, consult the detailed solver API docs. 
+`ITP()` is the recommended method for the scalar interval root-finding problems. It is particularly well-suited for cases where the function is smooth and well-behaved; and achieved superlinear convergence while retaining the optimal worst-case performance of the Bisection method. For more details, consult the detailed solver API docs.
 `Ridder` is a hybrid method that uses the value of function at the midpoint of the interval to perform an exponential interpolation to the root. This gives a fast convergence with a guaranteed convergence of at most twice the number of iterations as the bisection method.
 `Brent` is a combination of the bisection method, the secant method and inverse quadratic interpolation. At every iteration, Brent's method decides which method out of these three is likely to do best, and proceeds by doing a step according to that method. This gives a robust and fast method, which therefore enjoys considerable popularity.
 
diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md
index bf4ba3da4..6404d2226 100644
--- a/docs/src/tutorials/nonlinear.md
+++ b/docs/src/tutorials/nonlinear.md
@@ -47,6 +47,5 @@ p = 2.0
 probN = NonlinearProblem(f!, u0, p)
 
 linsolve = LinearSolve.KrylovJL_GMRES()
-sol = solve(probN, NewtonRaphson(), reltol = 1e-9, maxiters=4)
-sol = solve(probN, NewtonRaphson(; linsolve), reltol = 1e-9, maxiters=4)
+sol = solve(probN, NewtonRaphson(; linsolve), reltol = 1e-9)
 ```
diff --git a/test/basictests.jl b/test/basictests.jl
index 241faaf0f..a37b753f3 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -1,6 +1,7 @@
 using NonlinearSolve
 using StaticArrays
 using BenchmarkTools
+using LinearSolve
 using Test
 
 # --- NewtonRaphson tests ---

From 913231fcf9fce31aa35d2a92bb6df0ef3136a350 Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 22 Aug 2023 12:34:11 -0400
Subject: [PATCH 3/6] Reuse `jacobian_caches` & use `alg_autodiff` to choose
 correct `VecJac` operator

---
 src/NonlinearSolve.jl |  2 +-
 src/jacobian.jl       | 23 ++++++++++++++++++++++-
 src/levenberg.jl      | 20 --------------------
 src/raphson.jl        | 20 --------------------
 src/trustRegion.jl    | 20 --------------------
 5 files changed, 23 insertions(+), 62 deletions(-)

diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index 6a93a63e4..cae730bc3 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -34,10 +34,10 @@ function SciMLBase.__solve(prob::NonlinearProblem,
 end
 
 include("utils.jl")
-include("jacobian.jl")
 include("raphson.jl")
 include("trustRegion.jl")
 include("levenberg.jl")
+include("jacobian.jl")
 include("ad.jl")
 
 import PrecompileTools
diff --git a/src/jacobian.jl b/src/jacobian.jl
index 85c6479ad..8296069e0 100644
--- a/src/jacobian.jl
+++ b/src/jacobian.jl
@@ -92,7 +92,7 @@ function build_jac_and_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F
 
     J = if !linsolve_needs_jac
         # We don't need to construct the Jacobian
-        JacVec(uf, u; autodiff = AutoFiniteDiff())
+        JacVec(uf, u; autodiff = alg_autodiff(alg) ? AutoForwardDiff() : AutoFiniteDiff())
     else
         if f.jac_prototype === nothing
             ArrayInterface.undefmatrix(u)
@@ -104,6 +104,27 @@ function build_jac_and_jac_config(alg, f::F1, uf::F2, du1, u, tmp, du2) where {F
     return J, jac_config
 end
 
+# Build Jacobian Caches
+function jacobian_caches(alg::Union{NewtonRaphson, LevenbergMarquardt, TrustRegion}, f, u,
+    p, ::Val{true})
+    uf = JacobianWrapper(f, p)
+
+    du1 = zero(u)
+    du2 = zero(u)
+    tmp = zero(u)
+    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
+
+    linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
+    weight = similar(u)
+    recursivefill!(weight, true)
+
+    Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
+            nothing)..., weight)
+    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
+
+    uf, linsolve, J, du1, jac_config
+end
+
 function get_chunksize(jac_config::ForwardDiff.JacobianConfig{
     T,
     V,
diff --git a/src/levenberg.jl b/src/levenberg.jl
index f3a48f90f..db8955f4a 100644
--- a/src/levenberg.jl
+++ b/src/levenberg.jl
@@ -226,26 +226,6 @@ mutable struct LevenbergMarquardtCache{iip, fType, algType, uType, duType, resTy
     end
 end
 
-function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{true})
-    uf = JacobianWrapper(f, p)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
-
-    linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
-    weight = similar(u)
-    # Q: Setting this to false leads to residual = 0 in GMRES?
-    recursivefill!(weight, true)
-
-    Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
-            nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
-
-    uf, linsolve, J, du1, jac_config
-end
-
 function jacobian_caches(alg::LevenbergMarquardt, f, u, p, ::Val{false})
     JacobianWrapper(f, p), nothing, ArrayInterface.undefmatrix(u), nothing, nothing
 end
diff --git a/src/raphson.jl b/src/raphson.jl
index 572c68499..24e5799fd 100644
--- a/src/raphson.jl
+++ b/src/raphson.jl
@@ -104,26 +104,6 @@ mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, p
     end
 end
 
-function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true})
-    uf = JacobianWrapper(f, p)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
-
-    linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
-    weight = similar(u)
-    # Q: Setting this to false leads to residual = 0 in GMRES?
-    recursivefill!(weight, true)
-
-    Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
-            nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
-
-    uf, linsolve, J, du1, jac_config
-end
-
 function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{false})
     JacobianWrapper(f, p), nothing, nothing, nothing, nothing
 end
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
index 7c331cf97..6e867699c 100644
--- a/src/trustRegion.jl
+++ b/src/trustRegion.jl
@@ -278,26 +278,6 @@ mutable struct TrustRegionCache{iip, fType, algType, uType, resType, pType,
     end
 end
 
-function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{true})
-    uf = JacobianWrapper(f, p)
-
-    du1 = zero(u)
-    du2 = zero(u)
-    tmp = zero(u)
-    J, jac_config = build_jac_and_jac_config(alg, f, uf, du1, u, tmp, du2)
-
-    linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u)))
-    weight = similar(u)
-    # Q: Setting this to false leads to residual = 0 in GMRES?
-    recursivefill!(weight, true)
-
-    Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
-            nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
-
-    uf, linsolve, J, du1, jac_config
-end
-
 function jacobian_caches(alg::TrustRegion, f, u, p, ::Val{false})
     J = ArrayInterface.undefmatrix(u)
     JacobianWrapper(f, p), nothing, J, zero(u), nothing

From 6302183501e206c9520dea83ed70557869a207eb Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 22 Aug 2023 12:45:47 -0400
Subject: [PATCH 4/6] Add a test with Left Preconditioner

---
 test/basictests.jl | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/test/basictests.jl b/test/basictests.jl
index a37b753f3..1ab2fbaee 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -47,9 +47,9 @@ sol = benchmark_scalar(sf, csu0)
 # @test (@ballocated benchmark_mutable(ff, cu0)) < 200
 # @test (@ballocated benchmark_scalar(sf, csu0)) < 400
 
-function benchmark_inplace(f, u0, linsolve)
+function benchmark_inplace(f, u0, linsolve, precs)
     probN = NonlinearProblem{true}(f, u0)
-    solver = init(probN, NewtonRaphson(; linsolve), abstol = 1e-9)
+    solver = init(probN, NewtonRaphson(; linsolve, precs), abstol = 1e-9)
     sol = solve!(solver)
 end
 
@@ -58,8 +58,13 @@ function ffiip(du, u, p)
 end
 u0 = [1.0, 1.0]
 
-for linsolve in (nothing, KrylovJL_GMRES())
-    sol = benchmark_inplace(ffiip, u0, linsolve)
+precs = [
+    NonlinearSolve.DEFAULT_PRECS,
+    (args...) -> (Diagonal(rand!(similar(u0))), nothing)
+]
+
+for prec in precs, linsolve in (nothing, KrylovJL_GMRES())
+    sol = benchmark_inplace(ffiip, u0, linsolve, prec)
     @test sol.retcode === ReturnCode.Success
     @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
 end

From 0d989fb949629a660141c26ac9be749010f33dee Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 22 Aug 2023 13:06:46 -0400
Subject: [PATCH 5/6] Add random as test time dep

---
 Project.toml       | 3 ++-
 test/basictests.jl | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index f180547dc..3ae2b0bfc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -43,10 +43,11 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve"]
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random"]
diff --git a/test/basictests.jl b/test/basictests.jl
index 1ab2fbaee..76e1623e8 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -2,6 +2,7 @@ using NonlinearSolve
 using StaticArrays
 using BenchmarkTools
 using LinearSolve
+using Random
 using Test
 
 # --- NewtonRaphson tests ---

From c0067722275f3fad7b321ade9b51d35e6cbab68e Mon Sep 17 00:00:00 2001
From: Avik Pal <avikpal@mit.edu>
Date: Tue, 22 Aug 2023 13:15:09 -0400
Subject: [PATCH 6/6] Add LinearAlgebra as test time dep

---
 Project.toml       | 3 ++-
 test/basictests.jl | 1 +
 2 files changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 3ae2b0bfc..ed0b27f95 100644
--- a/Project.toml
+++ b/Project.toml
@@ -41,6 +41,7 @@ julia = "1.6"
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -50,4 +51,4 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random"]
+test = ["BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra"]
diff --git a/test/basictests.jl b/test/basictests.jl
index 76e1623e8..05a0152fa 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -3,6 +3,7 @@ using StaticArrays
 using BenchmarkTools
 using LinearSolve
 using Random
+using LinearAlgebra
 using Test
 
 # --- NewtonRaphson tests ---