SciML · ChrisRackauckas · Dec 25, 2023 · Dec 24, 2023 · Dec 24, 2023 · Dec 24, 2023
diff --git a/Project.toml b/Project.toml
@@ -1,7 +1,7 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 authors = ["SciML"]
-version = "3.3.0"
+version = "3.3.1"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"

diff --git a/docs/src/api/siamfanlequations.md b/docs/src/api/siamfanlequations.md
@@ -1,6 +1,7 @@
 # SIAMFANLEquations.jl
 
-This is an extension for importing solvers from [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) into the SciML
+This is an extension for importing solvers from
+[SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) into the SciML
 interface. Note that these solvers do not come by default, and thus one needs to install
 the package before using these solvers:
 

diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
@@ -14,8 +14,8 @@ solve(prob::SciMLBase.NonlinearProblem, args...; kwargs...)
 ## Iteration Controls
 
   - `maxiters::Int`: The maximum number of iterations to perform. Defaults to `1000`.
-  - `abstol::Number`: The absolute tolerance.
-  - `reltol::Number`: The relative tolerance.
+  - `abstol::Number`: The absolute tolerance. Defaults to `real(oneunit(T)) * (eps(real(one(T))))^(4 // 5)`.
+  - `reltol::Number`: The relative tolerance. Defaults to `real(oneunit(T)) * (eps(real(one(T))))^(2 // 5)`.
   - `termination_condition`: Termination Condition from DiffEqBase. Defaults to
     `AbsSafeBestTerminationMode()` for `NonlinearSolve.jl` and `AbsTerminateMode()` for
     `SimpleNonlinearSolve.jl`.

diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md
@@ -143,3 +143,10 @@ Newton-Krylov form. However, KINSOL is known to be less stable than some other
 implementations, as it has no line search or globalizer (trust region).
 
   - `KINSOL()`: The KINSOL method of the SUNDIALS C library
+
+### SIAMFANLEquations.jl
+
+SIAMFANLEquations.jl is a wrapper for the methods in the SIAMFANLEquations.jl library.
+
+  - `SIAMFANLEquationsJL()`: A wrapper for using the methods in
+    [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl)
diff --git a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl
@@ -5,7 +5,7 @@
 import FastLevenbergMarquardt as FastLM
 import FiniteDiff, ForwardDiff
 
-function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolve}
+@inline function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolve}
     if linsolve === :cholesky
         return FastLM.CholeskySolver(ArrayInterface.undefmatrix(x))
     elseif linsolve === :qr
@@ -33,14 +33,17 @@
 (f::InplaceFunction{false})(fx, x, p) = (fx .= f.f(x, p))
 
 function SciMLBase.__init(prob::NonlinearLeastSquaresProblem,
-        alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = 1e-8,
-        reltol = 1e-8, maxiters = 1000, kwargs...)
+        alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = nothing,
+        reltol = nothing, maxiters = 1000, kwargs...)
     iip = SciMLBase.isinplace(prob)
     u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
     fu = NonlinearSolve.evaluate_f(prob, u)
 
     f! = InplaceFunction{iip}(prob.f)
 
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
+    reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(u))
+
     if prob.f.jac === nothing
         use_forward_diff = if alg.autodiff === nothing
             ForwardDiff.can_dual(eltype(u))
@@ -95,9 +98,9 @@
     LM = FastLM.LMWorkspace(u, fu, J)
 
     return FastLevenbergMarquardtJLCache(f!, J!, prob, alg, LM, solver,
-        (; xtol = abstol, ftol = reltol, maxit = maxiters, alg.factor, alg.factoraccept,
-            alg.factorreject, alg.minscale, alg.maxscale, alg.factorupdate, alg.minfactor,
-            alg.maxfactor, kwargs...))
+        (; xtol = reltol, ftol = reltol, gtol = abstol, maxit = maxiters, alg.factor,
+            alg.factoraccept, alg.factorreject, alg.minscale, alg.maxscale,
+            alg.factorupdate, alg.minfactor, alg.maxfactor))
 end
 
 function SciMLBase.solve!(cache::FastLevenbergMarquardtJLCache)

diff --git a/ext/NonlinearSolveFixedPointAccelerationExt.jl b/ext/NonlinearSolveFixedPointAccelerationExt.jl
@@ -30,7 +30,7 @@
         f = (u) -> (prob.f(du, reshape(u, u_size), p); vec(du) .+ u)
     end
 
-    tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
     sol = fixed_point(f, NonlinearSolve._vec(u0); Algorithm = alg.algorithm,
         ConvergenceMetricThreshold = tol, MaxIter = maxiters, MaxM = alg.m,

diff --git a/ext/NonlinearSolveLeastSquaresOptimExt.jl b/ext/NonlinearSolveLeastSquaresOptimExt.jl
@@ -4,7 +4,7 @@
 import ConcreteStructs: @concrete
 import LeastSquaresOptim as LSO
 
-function _lso_solver(::LeastSquaresOptimJL{alg, linsolve}) where {alg, linsolve}
+@inline function _lso_solver(::LeastSquaresOptimJL{alg, linsolve}) where {alg, linsolve}
     ls = linsolve === :qr ? LSO.QR() :
          (linsolve === :cholesky ? LSO.Cholesky() :
           (linsolve === :lsmr ? LSO.LSMR() : nothing))
@@ -33,25 +33,28 @@
 (f::FunctionWrapper{false})(du, u) = (du .= f.f(u, f.p))
 
 function SciMLBase.__init(prob::NonlinearLeastSquaresProblem, alg::LeastSquaresOptimJL,
-        args...; alias_u0 = false, abstol = 1e-8, reltol = 1e-8, verbose = false,
-        maxiters = 1000, kwargs...)
+        args...; alias_u0 = false, abstol = nothing, show_trace::Val{ShT} = Val(false),
+        trace_level = TraceMinimal(), store_trace::Val{StT} = Val(false), maxiters = 1000,
+        reltol = nothing, kwargs...) where {ShT, StT}
     iip = SciMLBase.isinplace(prob)
-    u = alias_u0 ? prob.u0 : deepcopy(prob.u0)
+    u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
+
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
+    reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(u))
 
     f! = FunctionWrapper{iip}(prob.f, prob.p)
     g! = prob.f.jac === nothing ? nothing : FunctionWrapper{iip}(prob.f.jac, prob.p)
 
     resid_prototype = prob.f.resid_prototype === nothing ?
-                      (!iip ? prob.f(u, prob.p) : zeros(u)) :
-                      prob.f.resid_prototype
+                      (!iip ? prob.f(u, prob.p) : zeros(u)) : prob.f.resid_prototype
 
     lsoprob = LSO.LeastSquaresProblem(; x = u, f!, y = resid_prototype, g!,
         J = prob.f.jac_prototype, alg.autodiff, output_length = length(resid_prototype))
     allocated_prob = LSO.LeastSquaresProblemAllocated(lsoprob, _lso_solver(alg))
 
     return LeastSquaresOptimJLCache(prob, alg, allocated_prob,
-        (; x_tol = abstol, f_tol = reltol, iterations = maxiters, show_trace = verbose,
-            kwargs...))
+        (; x_tol = reltol, f_tol = abstol, g_tol = abstol, iterations = maxiters,
+            show_trace = ShT, store_trace = StT, show_every = trace_level.print_frequency))
 end
 
 function SciMLBase.solve!(cache::LeastSquaresOptimJLCache)

diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl
@@ -5,8 +5,9 @@
 
 function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
             NonlinearLeastSquaresProblem{uType, iip}}, alg::CMINPACK, args...;
-        abstol = nothing, maxiters = 100000, alias_u0::Bool = false,
-        termination_condition = nothing, kwargs...) where {uType, iip}
+        abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
+        show_trace::Val{ShT} = Val(false), store_trace::Val{StT} = Val(false),
+        termination_condition = nothing, kwargs...) where {uType, iip, ShT, StT}
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "CMINPACK does not support termination conditions!"
 
@@ -16,13 +17,12 @@
         u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
     end
 
-    T = eltype(u0)
     sizeu = size(prob.u0)
     p = prob.p
 
     # unwrapping alg params
-    show_trace = alg.show_trace
-    tracing = alg.tracing
+    show_trace = alg.show_trace || ShT
+    tracing = alg.tracing || StT
 
     if !iip && prob.u0 isa Number
         f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
@@ -44,7 +44,7 @@
     method = ifelse(alg.method === :auto,
         ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hybr), alg.method)
 
-    abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
 
     if SciMLBase.has_jac(prob.f)
         if !iip && prob.u0 isa Number
@@ -62,10 +62,10 @@
             end
         end
         original = MINPACK.fsolve(f!, g!, vec(u0), m; tol = abstol, show_trace, tracing,
-            method, iterations = maxiters, kwargs...)
+            method, iterations = maxiters)
     else
         original = MINPACK.fsolve(f!, vec(u0), m; tol = abstol, show_trace, tracing,
-            method, iterations = maxiters, kwargs...)
+            method, iterations = maxiters)
     end
 
     u = reshape(original.x, size(u))

diff --git a/ext/NonlinearSolveNLsolveExt.jl b/ext/NonlinearSolveNLsolveExt.jl
@@ -15,9 +15,7 @@
         u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
     end
 
-    T = eltype(u0)
     iip = isinplace(prob)
-
     sizeu = size(prob.u0)
     p = prob.p
 
@@ -70,7 +68,7 @@
         df = OnceDifferentiable(f!, vec(u0), vec(resid); autodiff)
     end
 
-    abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
     original = nlsolve(df, vec(u0); ftol = abstol, iterations = maxiters, method,
         store_trace, extended_trace, linesearch, linsolve, factor, autoscale, m, beta,

diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl
@@ -2,51 +2,58 @@
 
 using NonlinearSolve, SciMLBase
 using SIAMFANLEquations
-import ConcreteStructs: @concrete
 import UnPack: @unpack
 
-function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...; abstol = nothing,
-        reltol = nothing, alias_u0::Bool = false, maxiters = 1000, termination_condition = nothing, kwargs...)
-    @assert (termination_condition === nothing) || (termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
+@inline function __siam_fanl_equations_retcode_mapping(sol)
+    if sol.errcode == 0
+        return ReturnCode.Success
+    elseif sol.errcode == 10
+        return ReturnCode.MaxIters
+    elseif sol.errcode == 1
+        return ReturnCode.Failure
+    elseif sol.errcode == -1
+        return ReturnCode.Default
+    end
+end
+
+# pseudo transient continuation has a fixed cost per iteration, iteration statistics are
+# not interesting here.
+@inline function __siam_fanl_equations_stats_mapping(method, sol)
+    method === :pseudotransient && return nothing
+    return SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0,
+        sum(sol.stats.iarm))
+end
 
-    @unpack method, show_trace, delta, linsolve = alg
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...;
+        abstol = nothing, reltol = nothing, alias_u0::Bool = false, maxiters = 1000,
+        termination_condition = nothing, show_trace::Val{ShT} = Val(false),
+        kwargs...) where {ShT}
+    @assert (termination_condition ===
+             nothing)||(termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
+
+    @unpack method, delta, linsolve = alg
 
     iip = SciMLBase.isinplace(prob)
-    T = eltype(prob.u0)
 
-    atol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
-    rtol = reltol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : reltol
+    atol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(prob.u0))
+    rtol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(prob.u0))
 
     if prob.u0 isa Number
-        f! = if iip
-            function (u)
-                du = similar(u)
-                prob.f(du, u, prob.p)
-                return du
-            end
-        else
-            u -> prob.f(u, prob.p)
-        end
+        f = (u) -> prob.f(u, prob.p)
 
         if method == :newton
-            sol = nsolsc(f!, prob.u0; maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsolsc(f!, prob.u0; delta0 = delta, maxit = maxiters, atol = atol, rtol=rtol, printerr = show_trace)
+            sol = ptcsolsc(f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol,
+                printerr = ShT)
         elseif method == :secant
-            sol = secant(f!, prob.u0; maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
         end
 
-        if sol.errcode == 0
-            retcode = ReturnCode.Success
-        elseif sol.errcode == 10
-            retcode = ReturnCode.MaxIters
-        elseif sol.errcode == 1
-            retcode = ReturnCode.Failure
-        elseif sol.errcode == -1
-            retcode = ReturnCode.Default
-        end
-        stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
+        retcode = __siam_fanl_equations_retcode_mapping(sol)
+        stats = __siam_fanl_equations_stats_mapping(method, sol)
+        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode,
+            stats, original = sol)
     else
         u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
     end
@@ -71,67 +78,50 @@
     if linsolve !== nothing
         # Allocate ahead for Krylov basis
         JVS = linsolve == :gmres ? zeros(T, N, 3) : zeros(T, N)
-        # `linsolve` as a Symbol to keep unified interface with other EXTs, SIAMFANLEquations directly use String to choose between different linear solvers
+        # `linsolve` as a Symbol to keep unified interface with other EXTs,
+        # SIAMFANLEquations directly use String to choose between different linear solvers
         linsolve_alg = String(linsolve)
 
         if method == :newton
-            sol = nsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = nsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol,
+                rtol, printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
-        end
-
-        if sol.errcode == 0
-            retcode = ReturnCode.Success
-        elseif sol.errcode == 10
-            retcode = ReturnCode.MaxIters
-        elseif sol.errcode == 1
-            retcode = ReturnCode.Failure
-        elseif sol.errcode == -1
-            retcode = ReturnCode.Default
+            sol = ptcsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol,
+                rtol, printerr = ShT)
         end
-        stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
+
+        retcode = __siam_fanl_equations_retcode_mapping(sol)
+        stats = __siam_fanl_equations_stats_mapping(method, sol)
+        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode,
+            stats, original = sol)
     end
 
     # Allocate ahead for Jacobian
     FPS = zeros(T, N, N)
     if prob.f.jac === nothing
         # Use the built-in Jacobian machinery
         if method == :newton
-            sol = nsol(f!, u, FS, FPS;
-                        sham=1, atol = atol, rtol = rtol, maxit = maxiters,
-                        printerr = show_trace)
+            sol = nsol(f!, u, FS, FPS; sham = 1, atol, rtol, maxit = maxiters,
+                printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsol(f!, u, FS, FPS;
-                        atol = atol, rtol = rtol, maxit = maxiters,
-                        delta0 = delta, printerr = show_trace)
+            sol = ptcsol(f!, u, FS, FPS; atol, rtol, maxit = maxiters,
+                delta0 = delta, printerr = ShT)
         end
     else
         AJ!(J, u, x) = prob.f.jac(J, x, prob.p)
         if method == :newton
-            sol = nsol(f!, u, FS, FPS, AJ!;
-                        sham=1, atol = atol, rtol = rtol, maxit = maxiters,
-                        printerr = show_trace)
+            sol = nsol(f!, u, FS, FPS, AJ!; sham = 1, atol, rtol, maxit = maxiters,
+                printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsol(f!, u, FS, FPS, AJ!;
-                        atol = atol, rtol = rtol, maxit = maxiters,
-                        delta0 = delta, printerr = show_trace)
+            sol = ptcsol(f!, u, FS, FPS, AJ!; atol, rtol, maxit = maxiters,
+                delta0 = delta, printerr = ShT)
         end
     end
 
-    if sol.errcode == 0
-        retcode = ReturnCode.Success
-    elseif sol.errcode == 10
-        retcode = ReturnCode.MaxIters
-    elseif sol.errcode == 1
-        retcode = ReturnCode.Failure
-    elseif sol.errcode == -1
-        retcode = ReturnCode.Default
-    end
-
-    # pseudo transient continuation has a fixed cost per iteration, iteration statistics are not interesting here.
-    stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-    return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
+    retcode = __siam_fanl_equations_retcode_mapping(sol)
+    stats = __siam_fanl_equations_stats_mapping(method, sol)
+    return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats,
+        original = sol)
 end
 
-end
+end