preserve LM linearsolve structure

Author: avik-pal

.github/workflows/CI.yml

@@ -16,8 +16,7 @@ jobs:
     strategy:
       matrix:
         group:
-          - Core
-          - 23TestProblems
+          - All
         version:
           - '1'
     steps:

docs/src/solvers/NonlinearSystemSolvers.md

@@ -69,6 +69,8 @@ features, but have a bit of overhead on very small problems.
     Automatic Jacobian Resetting. This is a fast method but unstable for most problems!
   - `GeneralKlement()`: Generalization of Klement's Quasi-Newton Method with Line Search and
     Automatic Jacobian Resetting. This is a fast method but unstable for most problems!
+  - `LimitedMemoryBroyden()`: An advanced version of `LBroyden` which uses a limited memory
+    Broyden method. This is a fast method but unstable for most problems!
 
 ### SimpleNonlinearSolve.jl
 

src/jacobian.jl

@@ -50,8 +50,8 @@ jacobian!!(::Number, cache) = last(value_derivative(cache.uf, cache.u))
 
 # Build Jacobian Caches
 function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{iip};
-    linsolve_kwargs = (;),
-    linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ}
+    linsolve_kwargs = (;), lininit::Val{linsolve_init} = Val(true),
+    linsolve_with_JᵀJ::Val{needsJᵀJ} = Val(false)) where {iip, needsJᵀJ, linsolve_init}
     uf = JacobianWrapper{iip}(f, p)
 
     haslinsolve = hasfield(typeof(alg), :linsolve)
@@ -95,25 +95,28 @@ function jacobian_caches(alg::AbstractNonlinearSolveAlgorithm, f, u, p, ::Val{ii
         Jᵀfu = J' * _vec(fu)
     end
 
-    linprob = LinearProblem(needsJᵀJ ? __maybe_symmetric(JᵀJ) : J,
-        needsJᵀJ ? _vec(Jᵀfu) : _vec(fu); u0 = _vec(du))
-
-    if alg isa PseudoTransient
-        alpha = convert(eltype(u), alg.alpha_initial)
-        J_new = J - (1 / alpha) * I
-        linprob = LinearProblem(J_new, _vec(fu); u0 = _vec(du))
+    if linsolve_init
+        linprob_A = alg isa PseudoTransient ?
+                    (J - (1 / (convert(eltype(u), alg.alpha_initial))) * I) :
+                    (needsJᵀJ ? __maybe_symmetric(JᵀJ) : J)
+        linsolve = __setup_linsolve(linprob_A, needsJᵀJ ? Jᵀfu : fu, du, p, alg)
+    else
+        linsolve = nothing
     end
 
+    needsJᵀJ && return uf, linsolve, J, fu, jac_cache, du, JᵀJ, Jᵀfu
+    return uf, linsolve, J, fu, jac_cache, du
+end
+
+function __setup_linsolve(A, b, u, p, alg)
+    linprob = LinearProblem(A, _vec(b); u0 = _vec(u))
+
     weight = similar(u)
     recursivefill!(weight, true)
 
-    Pl, Pr = wrapprecs(alg.precs(needsJᵀJ ? __maybe_symmetric(JᵀJ) : J, nothing, u, p,
-            nothing, nothing, nothing, nothing, nothing)..., weight)
-    linsolve = init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr,
-        linsolve_kwargs...)
-
-    needsJᵀJ && return uf, linsolve, J, fu, jac_cache, du, JᵀJ, Jᵀfu
-    return uf, linsolve, J, fu, jac_cache, du
+    Pl, Pr = wrapprecs(alg.precs(A, nothing, u, p, nothing, nothing, nothing, nothing,
+            nothing)..., weight)
+    return init(linprob, alg.linsolve; alias_A = true, alias_b = true, Pl, Pr)
 end
 
 __get_nonsparse_ad(::AutoSparseForwardDiff) = AutoForwardDiff()

src/klement.jl

@@ -27,6 +27,10 @@ solves.
     linesearch
 end
 
+function set_linsolve(alg::GeneralKlement, linsolve)
+    return GeneralKlement(alg.max_resets, linsolve, alg.precs, alg.linesearch)
+end
+
 function GeneralKlement(; max_resets::Int = 5, linsolve = nothing,
     linesearch = LineSearch(), precs = DEFAULT_PRECS)
     linesearch = linesearch isa LineSearch ? linesearch : LineSearch(; method = linesearch)
@@ -60,7 +64,7 @@ end
 
 get_fu(cache::GeneralKlementCache) = cache.fu
 
-function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralKlement, args...;
+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
@@ -70,17 +74,13 @@ function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::GeneralKlemen
 
     if u isa Number
         linsolve = nothing
+        alg = alg_
     else
         # For General Julia Arrays default to LU Factorization
-        linsolve_alg = alg.linsolve === nothing && u isa Array ? LUFactorization() :
+        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...)
+        alg = set_linsolve(alg_, linsolve_alg)
+        linsolve = __setup_linsolve(J, _vec(fu), _vec(u), p, alg)
     end
 
     return GeneralKlementCache{iip}(f, alg, u, fu, zero(fu), _mutable_zero(u), p, linsolve,

src/levenberg.jl

@@ -213,10 +213,12 @@ function SciMLBase.__init(prob::Union{NonlinearProblem{uType, iip},
         mat_tmp = zero(JᵀJ)
         rhs_tmp = nothing
     else
-        mat_tmp = similar(JᵀJ, length(fu1) + length(u), length(u))
+        # Preserve Types
+        mat_tmp = vcat(J, DᵀD)
         fill!(mat_tmp, zero(eltype(u)))
-        rhs_tmp = similar(mat_tmp, length(fu1) + length(u))
+        rhs_tmp = vcat(fu1, u)
         fill!(rhs_tmp, zero(eltype(u)))
+        linsolve = __setup_linsolve(mat_tmp, rhs_tmp, u, p, alg)
     end
 
     return LevenbergMarquardtCache{iip, !_unwrap_val(linsolve_with_JᵀJ)}(f, alg, u, fu1,