diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index d3232532a..0d3ee1419 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -16,7 +16,8 @@ jobs:
     strategy:
       matrix:
         group:
-          - All
+          - Core
+          - 23TestProblems
         version:
           - '1'
     steps:
diff --git a/docs/pages.jl b/docs/pages.jl
index 8564261bd..f1b239bc9 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,8 +2,7 @@
 
 pages = ["index.md",
     "Getting Started with Nonlinear Rootfinding in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any[
-        "Code Optimization for Small Nonlinear Systems" => "tutorials/code_optimization.md",
+    "Tutorials" => Any["Code Optimization for Small Nonlinear Systems" => "tutorials/code_optimization.md",
         "Handling Large Ill-Conditioned and Sparse Systems" => "tutorials/large_systems.md",
         "Symbolic System Definition and Acceleration via ModelingToolkit" => "tutorials/modelingtoolkit.md",
         "tutorials/small_compile.md",
@@ -18,7 +17,8 @@ pages = ["index.md",
     "Solver Summaries and Recommendations" => Any["solvers/NonlinearSystemSolvers.md",
         "solvers/BracketingSolvers.md",
         "solvers/SteadyStateSolvers.md",
-        "solvers/NonlinearLeastSquaresSolvers.md"],
+        "solvers/NonlinearLeastSquaresSolvers.md",
+        "solvers/LineSearch.md"],
     "Detailed Solver APIs" => Any["api/nonlinearsolve.md",
         "api/simplenonlinearsolve.md",
         "api/minpack.md",
diff --git a/docs/src/api/nonlinearsolve.md b/docs/src/api/nonlinearsolve.md
index e4554f3e6..5da204f07 100644
--- a/docs/src/api/nonlinearsolve.md
+++ b/docs/src/api/nonlinearsolve.md
@@ -8,6 +8,9 @@ These are the native solvers of NonlinearSolve.jl.
 NewtonRaphson
 TrustRegion
 PseudoTransient
+DFSane
+GeneralBroyden
+GeneralKlement
 ```
 
 ## Polyalgorithms
diff --git a/docs/src/solvers/LineSearch.md b/docs/src/solvers/LineSearch.md
new file mode 100644
index 000000000..5d09301e2
--- /dev/null
+++ b/docs/src/solvers/LineSearch.md
@@ -0,0 +1,14 @@
+# [Line Search](@id linesearch)
+
+A convenience wrapper over `LineSearches.jl` and some native Line Search methods, powered
+internally with fast automatic differentiation.
+
+```@docs
+LineSearch
+```
+
+## Native Line Search Methods
+
+```@docs
+LiFukushimaLineSearch
+```
diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md
index de2d89a4e..7664bce10 100644
--- a/docs/src/solvers/NonlinearSystemSolvers.md
+++ b/docs/src/solvers/NonlinearSystemSolvers.md
@@ -65,6 +65,12 @@ features, but have a bit of overhead on very small problems.
   - `FastShortcutNonlinearPolyalg()`: The default method. A polyalgorithm that mixes fast methods
     with fallbacks to robust methods to allow for solving easy problems quickly without sacrificing
     robustnes on the hard problems.
+  - `GeneralBroyden()`: Generalization of Broyden's Quasi-Newton Method with Line Search and
+    Automatic Jacobian Resetting. This is a fast method but unstable when the condition number of
+    the Jacobian matrix is sufficiently large.
+  - `GeneralKlement()`: Generalization of Klement's Quasi-Newton Method with Line Search and
+    Automatic Jacobian Resetting. This is a fast method but unstable when the condition number of
+    the Jacobian matrix is sufficiently large.
 
 ### SimpleNonlinearSolve.jl
 
diff --git a/docs/src/tutorials/code_optimization.md b/docs/src/tutorials/code_optimization.md
index daefcb855..cd4089d5d 100644
--- a/docs/src/tutorials/code_optimization.md
+++ b/docs/src/tutorials/code_optimization.md
@@ -34,7 +34,7 @@ Take for example a prototypical small nonlinear solver code in its out-of-place
 ```@example small_opt
 using NonlinearSolve
 
-function f(u, p) 
+function f(u, p)
     u .* u .- p
 end
 u0 = [1.0, 1.0]
@@ -54,7 +54,7 @@ using BenchmarkTools
 Note that this way of writing the function is a shorthand for:
 
 ```@example small_opt
-function f(u, p) 
+function f(u, p)
     [u[1] * u[1] - p, u[2] * u[2] - p]
 end
 ```
@@ -69,25 +69,25 @@ In order to avoid this issue, we can use a non-allocating "in-place" approach. W
 this looks like:
 
 ```@example small_opt
-function f(du, u, p) 
+function f(du, u, p)
     du[1] = u[1] * u[1] - p
     du[2] = u[2] * u[2] - p
     nothing
 end
 
 prob = NonlinearProblem(f, u0, p)
-@btime  sol = solve(prob, NewtonRaphson())
+@btime sol = solve(prob, NewtonRaphson())
 ```
 
 Notice how much faster this already runs! We can make this code even simpler by using
 the `.=` in-place broadcasting.
 
 ```@example small_opt
-function f(du, u, p) 
+function f(du, u, p)
     du .= u .* u .- p
 end
 
-@btime  sol = solve(prob, NewtonRaphson())
+@btime sol = solve(prob, NewtonRaphson())
 ```
 
 ## Further Optimizations for Small Nonlinear Solves with Static Arrays and SimpleNonlinearSolve
@@ -140,7 +140,7 @@ want to use the out-of-place allocating form, but this time we want to output
 a static array. Doing it with broadcasting looks like:
 
 ```@example small_opt
-function f_SA(u, p) 
+function f_SA(u, p)
     u .* u .- p
 end
 u0 = SA[1.0, 1.0]
@@ -153,7 +153,7 @@ Note that only change here is that `u0` is made into a StaticArray! If we needed
 for a more complex nonlinear case, then we'd simply do the following:
 
 ```@example small_opt
-function f_SA(u, p) 
+function f_SA(u, p)
     SA[u[1] * u[1] - p, u[2] * u[2] - p]
 end
 
@@ -170,4 +170,4 @@ which are designed for these small-scale static examples. Let's now use `SimpleN
 @btime solve(prob, SimpleNewtonRaphson())
 ```
 
-And there we go, around 100ns from our starting point of almost 6μs!
\ No newline at end of file
+And there we go, around 100ns from our starting point of almost 6μs!
diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
index c1936932e..8e53d193c 100644
--- a/docs/src/tutorials/large_systems.md
+++ b/docs/src/tutorials/large_systems.md
@@ -1,8 +1,8 @@
 # [Efficiently Solving Large Sparse Ill-Conditioned Nonlinear Systems in Julia](@id large_systems)
 
-This tutorial is for getting into the extra features of using NonlinearSolve.jl. Solving ill-conditioned nonlinear systems 
-requires specializing the linear solver on properties of the Jacobian in order to cut down on the ``\mathcal{O}(n^3)`` 
-linear solve and the ``\mathcal{O}(n^2)`` back-solves. This tutorial is designed to explain the advanced usage of 
+This tutorial is for getting into the extra features of using NonlinearSolve.jl. Solving ill-conditioned nonlinear systems
+requires specializing the linear solver on properties of the Jacobian in order to cut down on the ``\mathcal{O}(n^3)``
+linear solve and the ``\mathcal{O}(n^2)`` back-solves. This tutorial is designed to explain the advanced usage of
 NonlinearSolve.jl by solving the steady state stiff Brusselator partial differential equation (BRUSS) using NonlinearSolve.jl.
 
 ## Definition of the Brusselator Equation
@@ -182,8 +182,8 @@ nothing # hide
 
 Notice that this acceleration does not require the definition of a sparsity
 pattern, and can thus be an easier way to scale for large problems. For more
-information on linear solver choices, see the 
-[linear solver documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#linear_nonlinear). 
+information on linear solver choices, see the
+[linear solver documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#linear_nonlinear).
 `linsolve` choices are any valid [LinearSolve.jl](https://linearsolve.sciml.ai/dev/) solver.
 
 !!! note
diff --git a/docs/src/tutorials/modelingtoolkit.md b/docs/src/tutorials/modelingtoolkit.md
index 7d26cb30c..3ee7d8e75 100644
--- a/docs/src/tutorials/modelingtoolkit.md
+++ b/docs/src/tutorials/modelingtoolkit.md
@@ -120,4 +120,4 @@ sol[u5]
 
 If you're interested in building models in a component or block based form, such as seen in systems like Simulink or Modelica,
 take a deeper look at [ModelingToolkit.jl's documentation](https://docs.sciml.ai/ModelingToolkit/stable/) which goes into
-detail on such features.
\ No newline at end of file
+detail on such features.
diff --git a/docs/src/tutorials/small_compile.md b/docs/src/tutorials/small_compile.md
index d7d2c2311..4f4c56b14 100644
--- a/docs/src/tutorials/small_compile.md
+++ b/docs/src/tutorials/small_compile.md
@@ -19,18 +19,18 @@ sol = solve(prob, SimpleNewtonRaphson())
 
 However, there are a few downsides to SimpleNonlinearSolve's `SimpleX` style algorithms to note:
 
-1. SimpleNonlinearSolve.jl's methods are not hooked into the LinearSolve.jl system, and thus do not have
-   the ability to specify linear solvers, use sparse matrices, preconditioners, and all of the other features
-   which are required to scale for very large systems of equations.
-2. SimpleNonlinearSolve.jl's methods have less robust error handling and termination conditions, and thus
-   these methods are missing some flexibility and give worse hints for debugging.
-3. SimpleNonlinearSolve.jl's methods are focused on out-of-place support. There is some in-place support,
-   but it's designed for simple smaller systems and so some optimizations are missing.
+ 1. SimpleNonlinearSolve.jl's methods are not hooked into the LinearSolve.jl system, and thus do not have
+    the ability to specify linear solvers, use sparse matrices, preconditioners, and all of the other features
+    which are required to scale for very large systems of equations.
+ 2. SimpleNonlinearSolve.jl's methods have less robust error handling and termination conditions, and thus
+    these methods are missing some flexibility and give worse hints for debugging.
+ 3. SimpleNonlinearSolve.jl's methods are focused on out-of-place support. There is some in-place support,
+    but it's designed for simple smaller systems and so some optimizations are missing.
 
 However, the major upsides of SimpleNonlinearSolve.jl are:
 
-1. The methods are optimized and non-allocating on StaticArrays
-2. The methods are minimal in compilation
+ 1. The methods are optimized and non-allocating on StaticArrays
+ 2. The methods are minimal in compilation
 
 As such, you can use the code as shown above to have very low startup with good methods, but for more scaling and debuggability
 we recommend the full NonlinearSolve.jl. But that said,
@@ -51,4 +51,4 @@ is not only sufficient but optimal.
 
 Julia has tools for building small binaries via static compilation with [StaticCompiler.jl](https://github.com/tshort/StaticCompiler.jl).
 However, these tools are currently limited to type-stable non-allocating functions. That said, SimpleNonlinearSolve.jl's solvers are
-precisely the subset of NonlinearSolve.jl which are compatible with static compilation.
\ No newline at end of file
+precisely the subset of NonlinearSolve.jl which are compatible with static compilation.
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index 08e5a7e79..2b26b3721 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -9,7 +9,7 @@ import ArrayInterface: restructure
 import ForwardDiff
 
 import ADTypes: AbstractFiniteDifferencesMode
-import ArrayInterface: undefmatrix, matrix_colors, parameterless_type, ismutable
+import ArrayInterface: undefmatrix, matrix_colors, parameterless_type, ismutable, issingular
 import ConcreteStructs: @concrete
 import EnumX: @enumx
 import ForwardDiff: Dual
@@ -26,6 +26,8 @@ import UnPack: @unpack
 const AbstractSparseADType = Union{ADTypes.AbstractSparseFiniteDifferences,
     ADTypes.AbstractSparseForwardMode, ADTypes.AbstractSparseReverseMode}
 
+abstract type AbstractNonlinearSolveLineSearchAlgorithm end
+
 abstract type AbstractNonlinearSolveAlgorithm <: AbstractNonlinearAlgorithm end
 abstract type AbstractNewtonAlgorithm{CJ, AD} <: AbstractNonlinearSolveAlgorithm end
 
@@ -50,10 +52,13 @@ function SciMLBase.solve!(cache::AbstractNonlinearSolveCache)
         cache.stats.nsteps += 1
     end
 
-    if cache.stats.nsteps == cache.maxiters
-        cache.retcode = ReturnCode.MaxIters
-    else
-        cache.retcode = ReturnCode.Success
+    # The solver might have set a different `retcode`
+    if cache.retcode == ReturnCode.Default
+        if cache.stats.nsteps == cache.maxiters
+            cache.retcode = ReturnCode.MaxIters
+        else
+            cache.retcode = ReturnCode.Success
+        end
     end
 
     return SciMLBase.build_solution(cache.prob, cache.alg, cache.u, get_fu(cache);
@@ -69,18 +74,22 @@ include("levenberg.jl")
 include("gaussnewton.jl")
 include("dfsane.jl")
 include("pseudotransient.jl")
+include("broyden.jl")
+include("klement.jl")
+include("lbroyden.jl")
 include("jacobian.jl")
 include("ad.jl")
 include("default.jl")
 
 import PrecompileTools
 
-@static if VERSION >= v"1.10"
+@static if VERSION ≥ v"1.10-"
     PrecompileTools.@compile_workload begin
         for T in (Float32, Float64)
             prob = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
 
-            precompile_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt())
+            precompile_algs = (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+                PseudoTransient(), GeneralBroyden(), GeneralKlement(), nothing)
 
             for alg in precompile_algs
                 solve(prob, alg, abstol = T(1e-2))
@@ -97,10 +106,11 @@ end
 
 export RadiusUpdateSchemes
 
-export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient
+export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
+    GeneralBroyden, GeneralKlement, LimitedMemoryBroyden
 export LeastSquaresOptimJL, FastLevenbergMarquardtJL
 export RobustMultiNewton, FastShortcutNonlinearPolyalg
 
-export LineSearch
+export LineSearch, LiFukushimaLineSearch
 
 end # module
diff --git a/src/broyden.jl b/src/broyden.jl
new file mode 100644
index 000000000..3232a2d9f
--- /dev/null
+++ b/src/broyden.jl
@@ -0,0 +1,166 @@
+# Sadly `Broyden` is taken up by SimpleNonlinearSolve.jl
+"""
+    GeneralBroyden(; max_resets = 3, linesearch = LineSearch(), reset_tolerance = nothing)
+
+An implementation of `Broyden` with reseting and line search.
+
+## Arguments
+
+  - `max_resets`: the maximum number of resets to perform. Defaults to `3`.
+  - `reset_tolerance`: the tolerance for the reset check. Defaults to
+    `sqrt(eps(eltype(u)))`.
+  - `linesearch`: the line search algorithm to use. Defaults to [`LineSearch()`](@ref),
+    which means that no line search is performed. Algorithms from `LineSearches.jl` can be
+    used here directly, and they will be converted to the correct `LineSearch`. It is
+    recommended to use [LiFukushimaLineSearchCache](@ref) -- a derivative free linesearch
+    specifically designed for Broyden's method.
+"""
+@concrete struct GeneralBroyden <: AbstractNewtonAlgorithm{false, Nothing}
+    max_resets::Int
+    reset_tolerance
+    linesearch
+end
+
+function GeneralBroyden(; max_resets = 3, linesearch = LineSearch(),
+    reset_tolerance = nothing)
+    linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
+    return GeneralBroyden(max_resets, reset_tolerance, linesearch)
+end
+
+@concrete mutable struct GeneralBroydenCache{iip} <: AbstractNonlinearSolveCache{iip}
+    f
+    alg
+    u
+    du
+    fu
+    fu2
+    dfu
+    p
+    J⁻¹
+    J⁻¹₂
+    J⁻¹df
+    force_stop::Bool
+    resets::Int
+    max_resets::Int
+    maxiters::Int
+    internalnorm
+    retcode::ReturnCode.T
+    abstol
+    reset_tolerance
+    reset_check
+    prob
+    stats::NLStats
+    lscache
+end
+
+get_fu(cache::GeneralBroydenCache) = cache.fu
+
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralBroyden, args...;
+    alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
+    kwargs...) where {uType, iip}
+    @unpack f, u0, p = prob
+    u = alias_u0 ? u0 : deepcopy(u0)
+    fu = evaluate_f(prob, u)
+    J⁻¹ = __init_identity_jacobian(u, fu)
+    reset_tolerance = alg.reset_tolerance === nothing ? sqrt(eps(eltype(u))) :
+                      alg.reset_tolerance
+    reset_check = x -> abs(x) ≤ reset_tolerance
+    return GeneralBroydenCache{iip}(f, alg, u, _mutable_zero(u), fu, zero(fu),
+        zero(fu), p, J⁻¹, zero(_vec(fu)'), _mutable_zero(u), false, 0, alg.max_resets,
+        maxiters, internalnorm, ReturnCode.Default, abstol, reset_tolerance,
+        reset_check, prob, NLStats(1, 0, 0, 0, 0),
+        init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)))
+end
+
+function perform_step!(cache::GeneralBroydenCache{true})
+    @unpack f, p, du, fu, fu2, dfu, u, J⁻¹, J⁻¹df, J⁻¹₂ = cache
+    T = eltype(u)
+
+    mul!(_vec(du), J⁻¹, -_vec(fu))
+    α = perform_linesearch!(cache.lscache, u, du)
+    axpy!(α, du, u)
+    f(fu2, u, p)
+
+    cache.internalnorm(fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the inverse jacobian
+    dfu .= fu2 .- fu
+
+    if all(cache.reset_check, du) || all(cache.reset_check, dfu)
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        fill!(J⁻¹, 0)
+        J⁻¹[diagind(J⁻¹)] .= T(1)
+        cache.resets += 1
+    else
+        mul!(_vec(J⁻¹df), J⁻¹, _vec(dfu))
+        mul!(J⁻¹₂, _vec(du)', J⁻¹)
+        du .= (du .- J⁻¹df) ./ (dot(du, J⁻¹df) .+ T(1e-5))
+        mul!(J⁻¹, _vec(du), J⁻¹₂, 1, 1)
+    end
+    fu .= fu2
+
+    return nothing
+end
+
+function perform_step!(cache::GeneralBroydenCache{false})
+    @unpack f, p = cache
+    T = eltype(cache.u)
+
+    cache.du = _restructure(cache.du, cache.J⁻¹ * -_vec(cache.fu))
+    α = perform_linesearch!(cache.lscache, cache.u, cache.du)
+    cache.u = cache.u .+ α * cache.du
+    cache.fu2 = f(cache.u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the inverse jacobian
+    cache.dfu = cache.fu2 .- cache.fu
+    if all(cache.reset_check, cache.du) || all(cache.reset_check, cache.dfu)
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        cache.J⁻¹ = __init_identity_jacobian(cache.u, cache.fu)
+        cache.resets += 1
+    else
+        cache.J⁻¹df = _restructure(cache.J⁻¹df, cache.J⁻¹ * _vec(cache.dfu))
+        cache.J⁻¹₂ = _vec(cache.du)' * cache.J⁻¹
+        cache.du = (cache.du .- cache.J⁻¹df) ./ (dot(cache.du, cache.J⁻¹df) .+ T(1e-5))
+        cache.J⁻¹ = cache.J⁻¹ .+ _vec(cache.du) * cache.J⁻¹₂
+    end
+    cache.fu = cache.fu2
+
+    return nothing
+end
+
+function SciMLBase.reinit!(cache::GeneralBroydenCache{iip}, u0 = cache.u; p = cache.p,
+    abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
+    cache.p = p
+    if iip
+        recursivecopy!(cache.u, u0)
+        cache.f(cache.fu, cache.u, p)
+    else
+        # don't have alias_u0 but cache.u is never mutated for OOP problems so it doesn't matter
+        cache.u = u0
+        cache.fu = cache.f(cache.u, p)
+    end
+    cache.abstol = abstol
+    cache.maxiters = maxiters
+    cache.stats.nf = 1
+    cache.stats.nsteps = 1
+    cache.resets = 0
+    cache.force_stop = false
+    cache.retcode = ReturnCode.Default
+    return cache
+end
diff --git a/src/default.jl b/src/default.jl
index 913f994fe..07effeecb 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -159,10 +159,8 @@ end
         ]
     else
         [
-            # FIXME: Broyden and Klement are type unstable
-            #        (upstream SimpleNonlinearSolve.jl issue)
-            !iip ? :(Klement()) : nothing, # Klement not yet implemented for IIP
-            !iip ? :(Broyden()) : nothing, # Broyden not yet implemented for IIP
+            :(GeneralBroyden()),
+            :(GeneralKlement()),
             :(NewtonRaphson(; linsolve, precs, adkwargs...)),
             :(NewtonRaphson(; linsolve, precs, linesearch = BackTracking(), adkwargs...)),
             :(TrustRegion(; linsolve, precs, adkwargs...)),
diff --git a/src/dfsane.jl b/src/dfsane.jl
index 13de5ff6a..f5bf69eca 100644
--- a/src/dfsane.jl
+++ b/src/dfsane.jl
@@ -43,7 +43,6 @@ See also the implementation in [SimpleNonlinearSolve.jl](https://github.com/SciM
 - `max_inner_iterations`: the maximum number of iterations allowed for the inner loop of the
   algorithm. Defaults to `1000`.
 """
-
 struct DFSane{T, F} <: AbstractNonlinearSolveAlgorithm
     σ_min::T
     σ_max::T
diff --git a/src/klement.jl b/src/klement.jl
new file mode 100644
index 000000000..e4d398445
--- /dev/null
+++ b/src/klement.jl
@@ -0,0 +1,217 @@
+"""
+    GeneralKlement(; max_resets = 5, linsolve = nothing,
+                     linesearch = LineSearch(), precs = DEFAULT_PRECS)
+
+An implementation of `Klement` with line search, preconditioning and customizable linear
+solves.
+
+## Keyword Arguments
+
+  - `max_resets`: the maximum number of resets to perform. Defaults to `5`.
+  - `linsolve`: the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) used for the
+    linear solves within the Newton method. Defaults to `nothing`, which means it uses the
+    LinearSolve.jl default algorithm choice. For more information on available algorithm
+    choices, see the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
+  - `precs`: the choice of preconditioners for the linear solver. Defaults to using no
+    preconditioners. For more information on specifying preconditioners for LinearSolve
+    algorithms, consult the
+    [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
+  - `linesearch`: the line search algorithm to use. Defaults to [`LineSearch()`](@ref),
+    which means that no line search is performed. Algorithms from `LineSearches.jl` can be
+    used here directly, and they will be converted to the correct `LineSearch`.
+"""
+@concrete struct GeneralKlement <: AbstractNewtonAlgorithm{false, Nothing}
+    max_resets::Int
+    linsolve
+    precs
+    linesearch
+end
+
+function GeneralKlement(; max_resets::Int = 5, linsolve = nothing,
+    linesearch = LineSearch(), precs = DEFAULT_PRECS)
+    linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
+    return GeneralKlement(max_resets, linsolve, precs, linesearch)
+end
+
+@concrete mutable struct GeneralKlementCache{iip} <: AbstractNonlinearSolveCache{iip}
+    f
+    alg
+    u
+    fu
+    fu2
+    du
+    p
+    linsolve
+    J
+    J_cache
+    J_cache2
+    Jᵀ²du
+    Jdu
+    resets
+    force_stop
+    maxiters::Int
+    internalnorm
+    retcode::ReturnCode.T
+    abstol
+    prob
+    stats::NLStats
+    lscache
+end
+
+get_fu(cache::GeneralKlementCache) = cache.fu
+
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralKlement, args...;
+    alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
+    linsolve_kwargs = (;), kwargs...) where {uType, iip}
+    @unpack f, u0, p = prob
+    u = alias_u0 ? u0 : deepcopy(u0)
+    fu = evaluate_f(prob, u)
+    J = __init_identity_jacobian(u, fu)
+
+    if u isa Number
+        linsolve = nothing
+    else
+        # For General Julia Arrays default to LU Factorization
+        linsolve_alg = alg.linsolve === nothing && u isa Array ? LUFactorization() :
+                       nothing
+        weight = similar(u)
+        recursivefill!(weight, true)
+        Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing,
+                nothing)..., weight)
+        linprob = LinearProblem(J, _vec(fu); u0 = _vec(fu))
+        linsolve = init(linprob, linsolve_alg; alias_A = true, alias_b = true, Pl, Pr,
+            linsolve_kwargs...)
+    end
+
+    return GeneralKlementCache{iip}(f, alg, u, fu, zero(fu), _mutable_zero(u), p, linsolve,
+        J, zero(J), zero(J), _vec(zero(fu)), _vec(zero(fu)), 0, false,
+        maxiters, internalnorm, ReturnCode.Default, abstol, prob, NLStats(1, 0, 0, 0, 0),
+        init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)))
+end
+
+function perform_step!(cache::GeneralKlementCache{true})
+    @unpack u, fu, f, p, alg, J, linsolve, du = cache
+    T = eltype(J)
+
+    singular, fact_done = _try_factorize_and_check_singular!(linsolve, J)
+
+    if singular
+        if cache.resets == alg.max_resets
+            cache.force_stop = true
+            cache.retcode = ReturnCode.Unstable
+            return nothing
+        end
+        fact_done = false
+        fill!(J, zero(T))
+        J[diagind(J)] .= T(1)
+        cache.resets += 1
+    end
+
+    # u = u - J \ fu
+    linres = dolinsolve(alg.precs, linsolve; A = ifelse(fact_done, nothing, J),
+        b = -_vec(fu), linu = _vec(du), p, reltol = cache.abstol)
+    cache.linsolve = linres.cache
+
+    # Line Search
+    α = perform_linesearch!(cache.lscache, u, du)
+    axpy!(α, du, u)
+    f(cache.fu2, u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+    cache.stats.nsolve += 1
+    cache.stats.nfactors += 1
+
+    cache.force_stop && return nothing
+
+    # Update the Jacobian
+    cache.J_cache .= cache.J' .^ 2
+    cache.Jdu .= _vec(du) .^ 2
+    mul!(cache.Jᵀ²du, cache.J_cache, cache.Jdu)
+    mul!(cache.Jdu, J, _vec(du))
+    cache.fu .= cache.fu2 .- cache.fu
+    cache.fu .= _restructure(cache.fu,
+        (_vec(cache.fu) .- cache.Jdu) ./ max.(cache.Jᵀ²du, eps(T)))
+    mul!(cache.J_cache, _vec(cache.fu), _vec(du)')
+    cache.J_cache .*= J
+    mul!(cache.J_cache2, cache.J_cache, J)
+    J .+= cache.J_cache2
+
+    cache.fu .= cache.fu2
+
+    return nothing
+end
+
+function perform_step!(cache::GeneralKlementCache{false})
+    @unpack fu, f, p, alg, J, linsolve = cache
+    T = eltype(J)
+
+    singular, fact_done = _try_factorize_and_check_singular!(linsolve, J)
+
+    if singular
+        if cache.resets == alg.max_resets
+            cache.force_stop = true
+            cache.retcode = ReturnCode.Unstable
+            return nothing
+        end
+        fact_done = false
+        cache.J = __init_identity_jacobian(cache.u, fu)
+        cache.resets += 1
+    end
+
+    # u = u - J \ fu
+    if linsolve === nothing
+        cache.du = -fu / cache.J
+    else
+        linres = dolinsolve(alg.precs, linsolve; A = ifelse(fact_done, nothing, J),
+            b = -_vec(fu), linu = _vec(cache.du), p, reltol = cache.abstol)
+        cache.linsolve = linres.cache
+    end
+
+    # Line Search
+    α = perform_linesearch!(cache.lscache, cache.u, cache.du)
+    cache.u = @. cache.u + α * cache.du  # `u` might not support mutation
+    cache.fu2 = f(cache.u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+    cache.stats.nsolve += 1
+    cache.stats.nfactors += 1
+
+    cache.force_stop && return nothing
+
+    # Update the Jacobian
+    cache.J_cache = cache.J' .^ 2
+    cache.Jdu = _vec(cache.du) .^ 2
+    cache.Jᵀ²du = cache.J_cache * cache.Jdu
+    cache.Jdu = J * _vec(cache.du)
+    cache.fu = cache.fu2 .- cache.fu
+    cache.fu = _restructure(cache.fu,
+        (_vec(cache.fu) .- cache.Jdu) ./ max.(cache.Jᵀ²du, eps(T)))
+    cache.J_cache = ((_vec(cache.fu) * _vec(cache.du)') .* J) * J
+    cache.J = J .+ cache.J_cache
+
+    cache.fu = cache.fu2
+
+    return nothing
+end
+
+function SciMLBase.reinit!(cache::GeneralKlementCache{iip}, u0 = cache.u; p = cache.p,
+    abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
+    cache.p = p
+    if iip
+        recursivecopy!(cache.u, u0)
+        cache.f(cache.fu, cache.u, p)
+    else
+        # don't have alias_u0 but cache.u is never mutated for OOP problems so it doesn't matter
+        cache.u = u0
+        cache.fu = cache.f(cache.u, p)
+    end
+    cache.abstol = abstol
+    cache.maxiters = maxiters
+    cache.stats.nf = 1
+    cache.stats.nsteps = 1
+    cache.force_stop = false
+    cache.retcode = ReturnCode.Default
+    return cache
+end
diff --git a/src/lbroyden.jl b/src/lbroyden.jl
new file mode 100644
index 000000000..458db8104
--- /dev/null
+++ b/src/lbroyden.jl
@@ -0,0 +1,260 @@
+"""
+    LimitedMemoryBroyden(; max_resets::Int = 3, linesearch = LineSearch(),
+        threshold::Int = 10, reset_tolerance = nothing)
+
+An implementation of `LimitedMemoryBroyden` with reseting and line search.
+
+## Arguments
+
+  - `max_resets`: the maximum number of resets to perform. Defaults to `3`.
+  - `reset_tolerance`: the tolerance for the reset check. Defaults to
+    `sqrt(eps(eltype(u)))`.
+  - `threshold`: the number of vectors to store in the low rank approximation. Defaults
+    to `10`.
+  - `linesearch`: the line search algorithm to use. Defaults to [`LineSearch()`](@ref),
+    which means that no line search is performed. Algorithms from `LineSearches.jl` can be
+    used here directly, and they will be converted to the correct `LineSearch`. It is
+    recommended to use [LiFukushimaLineSearchCache](@ref) -- a derivative free linesearch
+    specifically designed for Broyden's method.
+"""
+@concrete struct LimitedMemoryBroyden <: AbstractNewtonAlgorithm{false, Nothing}
+    max_resets::Int
+    threshold::Int
+    linesearch
+    reset_tolerance
+end
+
+function LimitedMemoryBroyden(; max_resets::Int = 3, linesearch = LineSearch(),
+    threshold::Int = 10, reset_tolerance = nothing)
+    linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
+    return LimitedMemoryBroyden(max_resets, threshold, linesearch, reset_tolerance)
+end
+
+@concrete mutable struct LimitedMemoryBroydenCache{iip} <: AbstractNonlinearSolveCache{iip}
+    f
+    alg
+    u
+    du
+    fu
+    fu2
+    dfu
+    p
+    U
+    Vᵀ
+    Ux
+    xᵀVᵀ
+    u_cache
+    vᵀ_cache
+    force_stop::Bool
+    resets::Int
+    iterations_since_reset::Int
+    max_resets::Int
+    maxiters::Int
+    internalnorm
+    retcode::ReturnCode.T
+    abstol
+    reset_tolerance
+    reset_check
+    prob
+    stats::NLStats
+    lscache
+end
+
+get_fu(cache::LimitedMemoryBroydenCache) = cache.fu
+
+function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::LimitedMemoryBroyden,
+    args...; alias_u0 = false, maxiters = 1000, abstol = 1e-6, internalnorm = DEFAULT_NORM,
+    kwargs...) where {uType, iip}
+    @unpack f, u0, p = prob
+    u = alias_u0 ? u0 : deepcopy(u0)
+    if u isa Number
+        # If u is a number then we simply use Broyden
+        return SciMLBase.__init(prob,
+            GeneralBroyden(; alg.max_resets, alg.reset_tolerance,
+                alg.linesearch), args...; alias_u0, maxiters, abstol, internalnorm, kwargs...)
+    end
+    fu = evaluate_f(prob, u)
+    threshold = min(alg.threshold, maxiters)
+    U, Vᵀ = __init_low_rank_jacobian(u, fu, threshold)
+    du = -fu
+    reset_tolerance = alg.reset_tolerance === nothing ? sqrt(eps(eltype(u))) :
+                      alg.reset_tolerance
+    reset_check = x -> abs(x) ≤ reset_tolerance
+    return LimitedMemoryBroydenCache{iip}(f, alg, u, du, fu, zero(fu),
+        zero(fu), p, U, Vᵀ, similar(u, threshold), similar(u, 1, threshold),
+        zero(u), zero(u), false, 0, 0, alg.max_resets, maxiters, internalnorm,
+        ReturnCode.Default, abstol, reset_tolerance, reset_check, prob,
+        NLStats(1, 0, 0, 0, 0),
+        init_linesearch_cache(alg.linesearch, f, u, p, fu, Val(iip)))
+end
+
+function perform_step!(cache::LimitedMemoryBroydenCache{true})
+    @unpack f, p, du, u = cache
+    T = eltype(u)
+
+    α = perform_linesearch!(cache.lscache, u, du)
+    axpy!(α, du, u)
+    f(cache.fu2, u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the Inverse Jacobian Approximation
+    cache.dfu .= cache.fu2 .- cache.fu
+
+    # Only try to reset if we have enough iterations since last reset
+    if cache.iterations_since_reset > size(cache.U, 1) &&
+       (all(cache.reset_check, du) || all(cache.reset_check, cache.dfu))
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        cache.iterations_since_reset = 0
+        cache.resets += 1
+        cache.du .= -cache.fu
+    else
+        idx = min(cache.iterations_since_reset, size(cache.U, 1))
+        U_part = selectdim(cache.U, 1, 1:idx)
+        Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx)
+
+        __lbroyden_matvec!(_vec(cache.vᵀ_cache), cache.Ux, U_part, Vᵀ_part, _vec(cache.du))
+        __lbroyden_rmatvec!(_vec(cache.u_cache), cache.xᵀVᵀ, U_part, Vᵀ_part,
+            _vec(cache.dfu))
+        cache.u_cache .= (du .- cache.u_cache) ./
+                         (dot(cache.vᵀ_cache, cache.dfu) .+ T(1e-5))
+
+        idx = mod1(cache.iterations_since_reset + 1, size(cache.U, 1))
+        selectdim(cache.U, 1, idx) .= _vec(cache.u_cache)
+        selectdim(cache.Vᵀ, 2, idx) .= _vec(cache.vᵀ_cache)
+
+        idx = min(cache.iterations_since_reset + 1, size(cache.U, 1))
+        U_part = selectdim(cache.U, 1, 1:idx)
+        Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx)
+        __lbroyden_matvec!(_vec(cache.du), cache.Ux, U_part, Vᵀ_part, _vec(cache.fu2))
+        cache.du .*= -1
+        cache.iterations_since_reset += 1
+    end
+
+    cache.fu .= cache.fu2
+
+    return nothing
+end
+
+function perform_step!(cache::LimitedMemoryBroydenCache{false})
+    @unpack f, p = cache
+    T = eltype(cache.u)
+
+    α = perform_linesearch!(cache.lscache, cache.u, cache.du)
+    cache.u = cache.u .+ α * cache.du
+    cache.fu2 = f(cache.u, p)
+
+    cache.internalnorm(cache.fu2) < cache.abstol && (cache.force_stop = true)
+    cache.stats.nf += 1
+
+    cache.force_stop && return nothing
+
+    # Update the Inverse Jacobian Approximation
+    cache.dfu .= cache.fu2 .- cache.fu
+
+    # Only try to reset if we have enough iterations since last reset
+    if cache.iterations_since_reset > size(cache.U, 1) &&
+       (all(cache.reset_check, cache.du) || all(cache.reset_check, cache.dfu))
+        if cache.resets ≥ cache.max_resets
+            cache.retcode = ReturnCode.Unstable
+            cache.force_stop = true
+            return nothing
+        end
+        cache.iterations_since_reset = 0
+        cache.resets += 1
+        cache.du = -cache.fu
+    else
+        idx = min(cache.iterations_since_reset, size(cache.U, 1))
+        U_part = selectdim(cache.U, 1, 1:idx)
+        Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx)
+
+        cache.vᵀ_cache = _restructure(cache.vᵀ_cache,
+            __lbroyden_matvec(U_part, Vᵀ_part, _vec(cache.du)))
+        cache.u_cache = _restructure(cache.u_cache,
+            __lbroyden_rmatvec(U_part, Vᵀ_part, _vec(cache.dfu)))
+        cache.u_cache = (cache.du .- cache.u_cache) ./
+                        (dot(cache.vᵀ_cache, cache.dfu) .+ T(1e-5))
+
+        idx = mod1(cache.iterations_since_reset + 1, size(cache.U, 1))
+        selectdim(cache.U, 1, idx) .= _vec(cache.u_cache)
+        selectdim(cache.Vᵀ, 2, idx) .= _vec(cache.vᵀ_cache)
+
+        idx = min(cache.iterations_since_reset + 1, size(cache.U, 1))
+        U_part = selectdim(cache.U, 1, 1:idx)
+        Vᵀ_part = selectdim(cache.Vᵀ, 2, 1:idx)
+        cache.du = _restructure(cache.du,
+            -__lbroyden_matvec(U_part, Vᵀ_part, _vec(cache.fu2)))
+        cache.iterations_since_reset += 1
+    end
+
+    cache.fu = cache.fu2
+
+    return nothing
+end
+
+function SciMLBase.reinit!(cache::LimitedMemoryBroydenCache{iip}, u0 = cache.u; p = cache.p,
+    abstol = cache.abstol, maxiters = cache.maxiters) where {iip}
+    cache.p = p
+    if iip
+        recursivecopy!(cache.u, u0)
+        cache.f(cache.fu, cache.u, p)
+    else
+        # don't have alias_u0 but cache.u is never mutated for OOP problems so it doesn't matter
+        cache.u = u0
+        cache.fu = cache.f(cache.u, p)
+    end
+    cache.abstol = abstol
+    cache.maxiters = maxiters
+    cache.stats.nf = 1
+    cache.stats.nsteps = 1
+    cache.resets = 0
+    cache.iterations_since_reset = 0
+    cache.force_stop = false
+    cache.retcode = ReturnCode.Default
+    return cache
+end
+
+@views function __lbroyden_matvec!(y::AbstractVector, Ux::AbstractVector,
+    U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector)
+    # Computes Vᵀ × U × x
+    η = size(U, 1)
+    if η == 0
+        y .= x
+        return nothing
+    end
+    mul!(Ux[1:η], U, x)
+    mul!(y, Vᵀ[:, 1:η], Ux[1:η])
+    return nothing
+end
+
+@views function __lbroyden_matvec(U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector)
+    # Computes Vᵀ × U × x
+    size(U, 1) == 0 && return x
+    return Vᵀ * (U * x)
+end
+
+@views function __lbroyden_rmatvec!(y::AbstractVector, xᵀVᵀ::AbstractMatrix,
+    U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector)
+    # Computes xᵀ × Vᵀ × U
+    η = size(U, 1)
+    if η == 0
+        y .= x
+        return nothing
+    end
+    mul!(xᵀVᵀ[:, 1:η], x', Vᵀ)
+    mul!(y', xᵀVᵀ[:, 1:η], U)
+    return nothing
+end
+
+@views function __lbroyden_rmatvec(U::AbstractMatrix, Vᵀ::AbstractMatrix, x::AbstractVector)
+    # Computes xᵀ × Vᵀ × U
+    size(U, 1) == 0 && return x
+    return (x' * Vᵀ) * U
+end
diff --git a/src/levenberg.jl b/src/levenberg.jl
index 016a12f31..fbfcc6a9a 100644
--- a/src/levenberg.jl
+++ b/src/levenberg.jl
@@ -301,7 +301,9 @@ function perform_step!(cache::LevenbergMarquardtCache{false})
     else
         linres = dolinsolve(alg.precs, linsolve;
             b = _mutable(_vec(J' *  #((2 / h) .* ((f(u .+ h .* v, p) .- fu1) ./ h .- J * v)))),
-                         _vec(((2 / h) .* ((_vec(f(u .+ h .* _restructure(u,v), p)) .- _vec(fu1)) ./ h .- J * _vec(v)))))),
+                              _vec(((2 / h) .*
+                                    ((_vec(f(u .+ h .* _restructure(u, v), p)) .-
+                                      _vec(fu1)) ./ h .- J * _vec(v)))))),
             linu = _vec(cache.a), p, reltol = cache.abstol)
         cache.linsolve = linres.cache
     end
@@ -311,7 +313,7 @@ function perform_step!(cache::LevenbergMarquardtCache{false})
     # Require acceptable steps to satisfy the following condition.
     norm_v = norm(v)
     if 2 * norm(cache.a) ≤ α_geodesic * norm_v
-        cache.δ = _restructure(cache.δ,_vec(v) .+ _vec(cache.a) ./ 2)
+        cache.δ = _restructure(cache.δ, _vec(v) .+ _vec(cache.a) ./ 2)
         @unpack δ, loss_old, norm_v_old, v_old, b_uphill = cache
         fu_new = f(u .+ δ, p)
         cache.stats.nf += 1
@@ -327,7 +329,7 @@ function perform_step!(cache::LevenbergMarquardtCache{false})
                 return nothing
             end
             cache.fu1 = fu_new
-            cache.v_old = _restructure(cache.v_old,v)
+            cache.v_old = _restructure(cache.v_old, v)
             cache.norm_v_old = norm_v
             cache.loss_old = loss
             cache.λ_factor = 1 / cache.damping_decrease_factor
diff --git a/src/linesearch.jl b/src/linesearch.jl
index 30861e14b..1f63d695c 100644
--- a/src/linesearch.jl
+++ b/src/linesearch.jl
@@ -26,7 +26,15 @@ function LineSearch(; method = Static(), autodiff = AutoFiniteDiff(), alpha = tr
     return LineSearch(method, autodiff, alpha)
 end
 
-@concrete mutable struct LineSearchCache
+@inline function init_linesearch_cache(ls::LineSearch, args...)
+    return init_linesearch_cache(ls.method, ls, args...)
+end
+
+# LineSearches.jl doesn't have a supertype so default to that
+init_linesearch_cache(_, ls, f, u, p, fu, iip) = LineSearchesJLCache(ls, f, u, p, fu, iip)
+
+# Wrapper over LineSearches.jl algorithms
+@concrete mutable struct LineSearchesJLCache
     f
     ϕ
     dϕ
@@ -35,11 +43,11 @@ end
     ls
 end
 
-function LineSearchCache(ls::LineSearch, f, u::Number, p, _, ::Val{false})
+function LineSearchesJLCache(ls::LineSearch, f, u::Number, p, _, ::Val{false})
     eval_f(u, du, α) = eval_f(u - α * du)
     eval_f(u) = f(u, p)
 
-    ls.method isa Static && return LineSearchCache(eval_f, nothing, nothing, nothing,
+    ls.method isa Static && return LineSearchesJLCache(eval_f, nothing, nothing, nothing,
         convert(typeof(u), ls.α), ls)
 
     g(u, fu) = last(value_derivative(Base.Fix2(f, p), u)) * fu
@@ -73,11 +81,11 @@ function LineSearchCache(ls::LineSearch, f, u::Number, p, _, ::Val{false})
         return ϕdϕ_internal
     end
 
-    return LineSearchCache(eval_f, ϕ, dϕ, ϕdϕ, convert(eltype(u), ls.α), ls)
+    return LineSearchesJLCache(eval_f, ϕ, dϕ, ϕdϕ, convert(eltype(u), ls.α), ls)
 end
 
-function LineSearchCache(ls::LineSearch, f, u, p, fu1, IIP::Val{iip}) where {iip}
-    fu = iip ? fu1 : nothing
+function LineSearchesJLCache(ls::LineSearch, f, u, p, fu1, IIP::Val{iip}) where {iip}
+    fu = iip ? deepcopy(fu1) : nothing
     u_ = _mutable_zero(u)
 
     function eval_f(u, du, α)
@@ -86,13 +94,14 @@ function LineSearchCache(ls::LineSearch, f, u, p, fu1, IIP::Val{iip}) where {iip
     end
     eval_f(u) = evaluate_f(f, u, p, IIP; fu)
 
-    ls.method isa Static && return LineSearchCache(eval_f, nothing, nothing, nothing,
+    ls.method isa Static && return LineSearchesJLCache(eval_f, nothing, nothing, nothing,
         convert(eltype(u), ls.α), ls)
 
     g₀ = _mutable_zero(u)
 
     autodiff = if iip && (ls.autodiff isa AutoZygote || ls.autodiff isa AutoSparseZygote)
-        @warn "Attempting to use Zygote.jl for linesearch on an in-place problem. Falling back to finite differencing."
+        @warn "Attempting to use Zygote.jl for linesearch on an in-place problem. Falling \
+               back to finite differencing."
         AutoFiniteDiff()
     else
         ls.autodiff
@@ -137,10 +146,10 @@ function LineSearchCache(ls::LineSearch, f, u, p, fu1, IIP::Val{iip}) where {iip
         return ϕdϕ_internal
     end
 
-    return LineSearchCache(eval_f, ϕ, dϕ, ϕdϕ, convert(eltype(u), ls.α), ls)
+    return LineSearchesJLCache(eval_f, ϕ, dϕ, ϕdϕ, convert(eltype(u), ls.α), ls)
 end
 
-function perform_linesearch!(cache::LineSearchCache, u, du)
+function perform_linesearch!(cache::LineSearchesJLCache, u, du)
     cache.ls.method isa Static && return cache.α
 
     ϕ = cache.ϕ(u, du)
@@ -154,3 +163,120 @@ function perform_linesearch!(cache::LineSearchCache, u, du)
 
     return first(cache.ls.method(ϕ, cache.dϕ(u, du), cache.ϕdϕ(u, du), cache.α, ϕ₀, dϕ₀))
 end
+
+"""
+    LiFukushimaLineSearch(; lambda_0 = 1.0, beta = 0.5, sigma_1 = 0.001,
+        eta = 0.1, nan_max_iter = 5, maxiters = 50)
+
+A derivative-free line search and global convergence of Broyden-like method for nonlinear
+equations by Dong-Hui Li & Masao Fukushima. For more details see
+https://doi.org/10.1080/10556780008805782
+"""
+struct LiFukushimaLineSearch{T} <: AbstractNonlinearSolveLineSearchAlgorithm
+    λ₀::T
+    β::T
+    σ₁::T
+    σ₂::T
+    η::T
+    ρ::T
+    nan_max_iter::Int
+    maxiters::Int
+end
+
+function LiFukushimaLineSearch(; lambda_0 = 1.0, beta = 0.1, sigma_1 = 0.001,
+    sigma_2 = 0.001, eta = 0.1, rho = 0.9, nan_max_iter = 5, maxiters = 50)
+    T = promote_type(typeof(lambda_0), typeof(beta), typeof(sigma_1), typeof(eta),
+        typeof(rho), typeof(sigma_2))
+    return LiFukushimaLineSearch{T}(lambda_0, beta, sigma_1, sigma_2, eta, rho,
+        nan_max_iter, maxiters)
+end
+
+@concrete mutable struct LiFukushimaLineSearchCache{iip}
+    f
+    p
+    u_cache
+    fu_cache
+    alg
+    α
+end
+
+function init_linesearch_cache(alg::LiFukushimaLineSearch, ls::LineSearch, f, _u, p, _fu,
+    ::Val{iip}) where {iip}
+    fu = iip ? deepcopy(_fu) : nothing
+    u = iip ? deepcopy(_u) : nothing
+    return LiFukushimaLineSearchCache{iip}(f, p, u, fu, alg, ls.α)
+end
+
+function perform_linesearch!(cache::LiFukushimaLineSearchCache{iip}, u, du) where {iip}
+    (; β, σ₁, σ₂, η, λ₀, ρ, nan_max_iter, maxiters) = cache.alg
+    λ₂ = λ₀
+    λ₁ = λ₂
+
+    if iip
+        cache.f(cache.fu_cache, u, cache.p)
+        fx_norm = norm(cache.fu_cache, 2)
+    else
+        fx_norm = norm(cache.f(u, cache.p), 2)
+    end
+
+    # Non-Blocking exit if the norm is NaN or Inf
+    !isfinite(fx_norm) && return cache.α
+
+    # Early Terminate based on Eq. 2.7
+    if iip
+        cache.u_cache .= u .+ du
+        cache.f(cache.fu_cache, cache.u_cache, cache.p)
+        fxλ_norm = norm(cache.fu_cache, 2)
+    else
+        fxλ_norm = norm(cache.f(u .+ du, cache.p), 2)
+    end
+
+    fxλ_norm ≤ ρ * fx_norm - σ₂ * norm(du, 2)^2 && return cache.α
+
+    if iip
+        cache.u_cache .= u .+ λ₂ .* du
+        cache.f(cache.fu_cache, cache.u_cache, cache.p)
+        fxλp_norm = norm(cache.fu_cache, 2)
+    else
+        fxλp_norm = norm(cache.f(u .+ λ₂ .* du, cache.p), 2)
+    end
+
+    if !isfinite(fxλp_norm)
+        # Backtrack a finite number of steps
+        nan_converged = false
+        for _ in 1:nan_max_iter
+            λ₁, λ₂ = λ₂, β * λ₂
+
+            if iip
+                cache.u_cache .= u .+ λ₂ .* du
+                cache.f(cache.fu_cache, cache.u_cache, cache.p)
+                fxλp_norm = norm(cache.fu_cache, 2)
+            else
+                fxλp_norm = norm(cache.f(u .+ λ₂ .* du, cache.p), 2)
+            end
+
+            nan_converged = isfinite(fxλp_norm)
+            nan_converged && break
+        end
+
+        # Non-Blocking exit if the norm is still NaN or Inf
+        !nan_converged && return cache.α
+    end
+
+    for _ in 1:maxiters
+        if iip
+            cache.u_cache .= u .+ λ₂ .* du
+            cache.f(cache.fu_cache, cache.u_cache, cache.p)
+            fxλp_norm = norm(cache.fu_cache, 2)
+        else
+            fxλp_norm = norm(cache.f(u .+ λ₂ .* du, cache.p), 2)
+        end
+
+        converged = fxλp_norm ≤ (1 + η) * fx_norm - σ₁ * λ₂^2 * norm(du, 2)^2
+
+        converged && break
+        λ₁, λ₂ = λ₂, β * λ₂
+    end
+
+    return λ₂
+end
diff --git a/src/pseudotransient.jl b/src/pseudotransient.jl
index 8e70f2350..a041c258c 100644
--- a/src/pseudotransient.jl
+++ b/src/pseudotransient.jl
@@ -29,10 +29,6 @@ SIAM Journal on Scientific Computing,25, 553-569.](https://doi.org/10.1137/S1064
     [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/).
   - `alpha_initial` : the initial pseudo time step. it defaults to 1e-3. If it is small, 
     you are going to need more iterations to converge but it can be more stable.
-
-
-
-
 """
 @concrete struct PseudoTransient{CJ, AD} <: AbstractNewtonAlgorithm{CJ, AD}
     ad::AD
@@ -92,19 +88,14 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::PseudoTransi
     else
         fu1 = _mutable(f(u, p))
     end
-    uf, linsolve, J, fu2, jac_cache, du = jacobian_caches(alg,
-        f,
-        u,
-        p,
-        Val(iip);
+    uf, linsolve, J, fu2, jac_cache, du = jacobian_caches(alg, f, u, p, Val(iip);
         linsolve_kwargs)
     alpha = convert(eltype(u), alg.alpha_initial)
     res_norm = internalnorm(fu1)
 
     return PseudoTransientCache{iip}(f, alg, u, fu1, fu2, du, p, alpha, res_norm, uf,
-        linsolve, J,
-        jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob,
-        NLStats(1, 0, 0, 0, 0))
+        linsolve, J, jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol,
+        prob, NLStats(1, 0, 0, 0, 0))
 end
 
 function perform_step!(cache::PseudoTransientCache{true})
diff --git a/src/raphson.jl b/src/raphson.jl
index 642c7ee36..db9e9e322 100644
--- a/src/raphson.jl
+++ b/src/raphson.jl
@@ -82,7 +82,8 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg_::NewtonRaphso
 
     return NewtonRaphsonCache{iip}(f, alg, u, fu1, fu2, du, p, uf, linsolve, J,
         jac_cache, false, maxiters, internalnorm, ReturnCode.Default, abstol, prob,
-        NLStats(1, 0, 0, 0, 0), LineSearchCache(alg.linesearch, f, u, p, fu1, Val(iip)))
+        NLStats(1, 0, 0, 0, 0),
+        init_linesearch_cache(alg.linesearch, f, u, p, fu1, Val(iip)))
 end
 
 function perform_step!(cache::NewtonRaphsonCache{true})
@@ -96,7 +97,7 @@ function perform_step!(cache::NewtonRaphsonCache{true})
 
     # Line Search
     α = perform_linesearch!(cache.lscache, u, du)
-    @. u = u - α * du
+    axpy!(-α, du, u)
     f(cache.fu1, u, p)
 
     cache.internalnorm(fu1) < cache.abstol && (cache.force_stop = true)
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
index f97b63637..8e7118cc6 100644
--- a/src/trustRegion.jl
+++ b/src/trustRegion.jl
@@ -81,7 +81,6 @@ states as `RadiusUpdateSchemes.T`. Simply put the desired scheme as follows:
 end
 
 """
-```julia
     TrustRegion(; concrete_jac = nothing, linsolve = nothing, precs = DEFAULT_PRECS,
         radius_update_scheme::RadiusUpdateSchemes.T = RadiusUpdateSchemes.Simple,
         max_trust_radius::Real = 0 // 1, initial_trust_radius::Real = 0 // 1,
@@ -419,7 +418,8 @@ function trust_region_step!(cache::TrustRegionCache)
     cache.loss_new = get_loss(fu_new)
 
     # Compute the ratio of the actual reduction to the predicted reduction.
-    cache.r = -(loss - cache.loss_new) / (dot(_vec(du), _vec(g)) + dot(_vec(du), H, _vec(du)) / 2)
+    cache.r = -(loss - cache.loss_new) /
+              (dot(_vec(du), _vec(g)) + dot(_vec(du), H, _vec(du)) / 2)
     @unpack r = cache
 
     if radius_update_scheme === RadiusUpdateSchemes.Simple
diff --git a/src/utils.jl b/src/utils.jl
index a26d1862a..c99b05bf4 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -35,8 +35,8 @@ function default_adargs_to_adtype(; chunk_size = missing, autodiff = nothing,
     if chunk_size !== missing || standardtag !== missing || diff_type !== missing ||
        autodiff !== missing
         Base.depwarn("`chunk_size`, `standardtag`, `diff_type`, \
-            `autodiff::Union{Val, Bool}` kwargs have been deprecated and will be removed in \
-             v3. Update your code to directly specify autodiff=<ADTypes>",
+            `autodiff::Union{Val, Bool}` kwargs have been deprecated and will be removed \
+             in v3. Update your code to directly specify autodiff=<ADTypes>",
             :default_adargs_to_adtype)
     end
     chunk_size === missing && (chunk_size = Val{0}())
@@ -74,8 +74,8 @@ end
 @inline _vec(v::Number) = v
 @inline _vec(v::AbstractVector) = v
 
-@inline _restructure(y,x) = restructure(y,x)
-@inline _restructure(y::Number,x::Number) = x
+@inline _restructure(y, x) = restructure(y, x)
+@inline _restructure(y::Number, x::Number) = x
 
 DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing
 
@@ -211,3 +211,50 @@ function __get_concrete_algorithm(alg, prob)
     end
     return set_ad(alg, ad)
 end
+
+__init_identity_jacobian(u::Number, _) = u
+function __init_identity_jacobian(u, fu)
+    return convert(parameterless_type(_mutable(u)),
+        Matrix{eltype(u)}(I, length(fu), length(u)))
+end
+function __init_identity_jacobian(u::StaticArray, fu)
+    return convert(MArray{Tuple{length(fu), length(u)}},
+        Matrix{eltype(u)}(I, length(fu), length(u)))
+end
+
+function __init_low_rank_jacobian(u::StaticArray, fu, threshold::Int)
+    Vᵀ = convert(MArray{Tuple{length(u), threshold}},
+        zeros(eltype(u), length(u), threshold))
+    U = convert(MArray{Tuple{threshold, length(u)}}, zeros(eltype(u), threshold, length(u)))
+    return U, Vᵀ
+end
+function __init_low_rank_jacobian(u, fu, threshold::Int)
+    Vᵀ = convert(parameterless_type(_mutable(u)), zeros(eltype(u), length(u), threshold))
+    U = convert(parameterless_type(_mutable(u)), zeros(eltype(u), threshold, length(u)))
+    return U, Vᵀ
+end
+
+# Check Singular Matrix
+_issingular(x::Number) = iszero(x)
+@generated function _issingular(x::T) where {T}
+    hasmethod(issingular, Tuple{T}) && return :(issingular(x))
+    return :(__issingular(x))
+end
+__issingular(x::AbstractMatrix{T}) where {T} = cond(x) > inv(sqrt(eps(T)))
+__issingular(x) = false ## If SciMLOperator and such
+
+# If factorization is LU then perform that and update the linsolve cache
+# else check if the matrix is singular
+function _try_factorize_and_check_singular!(linsolve, X)
+    if linsolve.cacheval isa LU
+        # LU Factorization was used
+        linsolve.A = X
+        linsolve.cacheval = LinearSolve.do_factorization(linsolve.alg, X, linsolve.b,
+            linsolve.u)
+        linsolve.isfresh = false
+
+        return !issuccess(linsolve.cacheval), true
+    end
+    return _issingular(X), false
+end
+_try_factorize_and_check_singular!(::Nothing, x) = _issingular(x), false
diff --git a/test/23_test_problems.jl b/test/23_test_problems.jl
index 5e0df8352..2bdbecffb 100644
--- a/test/23_test_problems.jl
+++ b/test/23_test_problems.jl
@@ -8,19 +8,27 @@ function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4)
         x = dict["start"]
         res = similar(x)
         nlprob = NonlinearProblem(problem, x)
-        @testset "$(dict["title"])" begin
+        @testset "$idx: $(dict["title"])" begin
             for alg in alg_ops
-                sol = solve(nlprob, alg, abstol = 1e-18, reltol = 1e-18)
-                problem(res, sol.u, nothing)
-                broken = idx in broken_tests[alg] ? true : false
-                @test norm(res)≤ϵ broken=broken
+                try
+                    sol = solve(nlprob, alg, abstol = 1e-18, reltol = 1e-18)
+                    problem(res, sol.u, nothing)
+                    broken = idx in broken_tests[alg] ? true : false
+                    @test norm(res)≤ϵ broken=broken
+                catch
+                    broken = idx in broken_tests[alg] ? true : false
+                    if broken
+                        @test false broken=true
+                    else
+                        @test 1 == 2
+                    end
+                end
             end
         end
     end
 end
 
-# NewtonRaphson
-@testset "NewtonRaphson test problem library" begin
+@testset "NewtonRaphson 23 Test Problems" begin
     alg_ops = (NewtonRaphson(),)
 
     # dictionary with indices of test problems where method does not converge to small residual
@@ -30,7 +38,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "TrustRegion test problem library" begin
+@testset "TrustRegion 23 Test Problems" begin
     alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple),
         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan),
         TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei),
@@ -50,7 +58,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-@testset "TrustRegion test problem library" begin
+@testset "LevenbergMarquardt 23 Test Problems" begin
     alg_ops = (LevenbergMarquardt(; linsolve = NormalCholeskyFactorization()),
         LevenbergMarquardt(; α_geodesic = 0.1, linsolve = NormalCholeskyFactorization()))
 
@@ -62,8 +70,7 @@ end
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
 
-# DFSane
-@testset "DFSane test problem library" begin
+@testset "DFSane 23 Test Problems" begin
     alg_ops = (DFSane(),)
 
     broken_tests = Dict(alg => Int[] for alg in alg_ops)
@@ -71,3 +78,31 @@ end
 
     test_on_library(problems, dicts, alg_ops, broken_tests)
 end
+
+@testset "GeneralBroyden 23 Test Problems" begin
+    alg_ops = (GeneralBroyden(),
+        GeneralBroyden(; linesearch = LiFukushimaLineSearch(; beta = 0.1)),
+        GeneralBroyden(; linesearch = BackTracking()))
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [1, 3, 4, 5, 6, 11, 12, 13, 14, 21]
+    broken_tests[alg_ops[2]] = [1, 2, 3, 4, 5, 6, 9, 11, 13, 22]
+    broken_tests[alg_ops[3]] = [1, 2, 4, 5, 6, 11, 12, 13, 14, 21]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+@testset "GeneralKlement 23 Test Problems" begin
+    alg_ops = (GeneralKlement(),
+        GeneralKlement(; linesearch = BackTracking()),
+        GeneralKlement(; linesearch = HagerZhang()))
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 7, 11, 13, 22]
+    broken_tests[alg_ops[2]] = [1, 2, 4, 5, 6, 7, 11, 12, 13, 22]
+    broken_tests[alg_ops[3]] = [1, 2, 4, 5, 6, 8, 11, 12, 13, 22]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+# NOTE: Not adding LimitedMemoryBroyden here since it fails on most of the preoblems
diff --git a/test/basictests.jl b/test/basictests.jl
index f20417785..681d506de 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -664,3 +664,272 @@ end
         @test all(abs.(newton_fails(sol.u, p)) .< 1e-10)
     end
 end
+
+# --- GeneralBroyden tests ---
+
+@testset "GeneralBroyden" begin
+    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, GeneralBroyden(; linesearch), abstol = 1e-9)
+    end
+
+    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{true}(f, u0, p)
+        return solve(prob, GeneralBroyden(; linesearch), abstol = 1e-9)
+    end
+
+    @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
+            StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(),
+            LiFukushimaLineSearch()),
+        ad in (AutoFiniteDiff(), AutoZygote())
+
+        linesearch = LineSearch(; method = lsmethod, autodiff = ad)
+        u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0),
+                GeneralBroyden(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) < 200
+        end
+
+        @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+            ad isa AutoZygote && continue
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
+                GeneralBroyden(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) ≤ 64
+        end
+    end
+
+    @testset "[OOP] [Immutable AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+                res_true = sqrt(p)
+                all(res.u .≈ res_true)
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                @SVector[1.0, 1.0], p).u[end], p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
+
+    @testset "[OOP] [Scalar AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+                res_true = sqrt(p)
+                res.u ≈ res_true
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u,
+                p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
+
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
+
+    # Iterator interface
+    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
+        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
+        cache = init(probN, GeneralBroyden(); maxiters = 100, abstol = 1e-10)
+        sols = zeros(length(p_range))
+        for (i, p) in enumerate(p_range)
+            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
+            sol = solve!(cache)
+            sols[i] = iip ? sol.u[1] : sol.u
+        end
+        return sols
+    end
+    p = range(0.01, 2, length = 200)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+end
+
+# --- GeneralKlement tests ---
+
+@testset "GeneralKlement" begin
+    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, GeneralKlement(; linesearch), abstol = 1e-9)
+    end
+
+    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{true}(f, u0, p)
+        return solve(prob, GeneralKlement(; linesearch), abstol = 1e-9)
+    end
+
+    @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
+            StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente()),
+        ad in (AutoFiniteDiff(), AutoZygote())
+
+        linesearch = LineSearch(; method = lsmethod, autodiff = ad)
+        u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0),
+                GeneralKlement(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) < 200
+        end
+
+        @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+            ad isa AutoZygote && continue
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
+                GeneralKlement(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) ≤ 64
+        end
+    end
+
+    @testset "[OOP] [Immutable AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+                res_true = sqrt(p)
+                all(res.u .≈ res_true)
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                @SVector[1.0, 1.0], p).u[end], p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
+
+    @testset "[OOP] [Scalar AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+                res_true = sqrt(p)
+                res.u ≈ res_true
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u,
+                p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
+
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
+
+    # Iterator interface
+    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
+        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
+        cache = init(probN, GeneralKlement(); maxiters = 100, abstol = 1e-10)
+        sols = zeros(length(p_range))
+        for (i, p) in enumerate(p_range)
+            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
+            sol = solve!(cache)
+            sols[i] = iip ? sol.u[1] : sol.u
+        end
+        return sols
+    end
+    p = range(0.01, 2, length = 200)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false)) ≈ sqrt.(p)
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true)) ≈ sqrt.(p)
+end
+
+# --- LimitedMemoryBroyden tests ---
+
+@testset "LimitedMemoryBroyden" begin
+    function benchmark_nlsolve_oop(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, LimitedMemoryBroyden(; linesearch), abstol = 1e-9)
+    end
+
+    function benchmark_nlsolve_iip(f, u0, p = 2.0; linesearch = LineSearch())
+        prob = NonlinearProblem{true}(f, u0, p)
+        return solve(prob, LimitedMemoryBroyden(; linesearch), abstol = 1e-9)
+    end
+
+    @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (Static(),
+            StrongWolfe(), BackTracking(), HagerZhang(), MoreThuente(),
+            LiFukushimaLineSearch()),
+        ad in (AutoFiniteDiff(), AutoZygote())
+
+        linesearch = LineSearch(; method = lsmethod, autodiff = ad)
+        u0s = ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in u0s
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{false}(quadratic_f, u0, 2.0),
+                LimitedMemoryBroyden(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) < 200
+        end
+
+        @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+            ad isa AutoZygote && continue
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; linesearch)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+
+            cache = init(NonlinearProblem{true}(quadratic_f!, u0, 2.0),
+                LimitedMemoryBroyden(; linesearch), abstol = 1e-9)
+            @test (@ballocated solve!($cache)) ≤ 64
+        end
+    end
+
+    @testset "[OOP] [Immutable AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+                res_true = sqrt(p)
+                all(((x, y),) -> isapprox(x, y; atol = 1e-3), zip(res.u, res_true))
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                @SVector[1.0, 1.0], p).u[end], p)≈1 / (2 * sqrt(p)) atol=1e-3
+        end
+    end
+
+    @testset "[OOP] [Scalar AD]" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+                res_true = sqrt(p)
+                res.u ≈ res_true
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u,
+                p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
+
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
+
+    # Iterator interface
+    function nlprob_iterator_interface(f, p_range, ::Val{iip}) where {iip}
+        probN = NonlinearProblem{iip}(f, iip ? [0.5] : 0.5, p_range[begin])
+        cache = init(probN, LimitedMemoryBroyden(); maxiters = 100, abstol = 1e-10)
+        sols = zeros(length(p_range))
+        for (i, p) in enumerate(p_range)
+            reinit!(cache, iip ? [cache.u[1]] : cache.u; p = p)
+            sol = solve!(cache)
+            sols[i] = iip ? sol.u[1] : sol.u
+        end
+        return sols
+    end
+    p = range(0.01, 2, length = 200)
+    @test nlprob_iterator_interface(quadratic_f, p, Val(false))≈sqrt.(p) atol=1e-2
+    @test nlprob_iterator_interface(quadratic_f!, p, Val(true))≈sqrt.(p) atol=1e-2
+end
diff --git a/test/matrix_resizing.jl b/test/matrix_resizing.jl
index f8ad423e5..1d9462fa1 100644
--- a/test/matrix_resizing.jl
+++ b/test/matrix_resizing.jl
@@ -1,21 +1,24 @@
 using NonlinearSolve, Test
 
 ff(u, p) = u .* u .- p
-u0 = rand(2,2)
+u0 = rand(2, 2)
 p = 2.0
 vecprob = NonlinearProblem(ff, vec(u0), p)
 prob = NonlinearProblem(ff, u0, p)
 
-for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg())
+for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
+    RobustMultiNewton(), FastShortcutNonlinearPolyalg(), GeneralBroyden(), GeneralKlement(),
+    LimitedMemoryBroyden())
     @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
 end
 
 fiip(du, u, p) = (du .= u .* u .- p)
-u0 = rand(2,2)
+u0 = rand(2, 2)
 p = 2.0
 vecprob = NonlinearProblem(fiip, vec(u0), p)
 prob = NonlinearProblem(fiip, u0, p)
 
-for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg())
+for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
+    RobustMultiNewton(), FastShortcutNonlinearPolyalg(), GeneralBroyden(), GeneralKlement())
     @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
-end
\ No newline at end of file
+end