Standardize parts of SIAM FANL Equations

Author: avik-pal

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:
 

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)

ext/NonlinearSolveNLsolveExt.jl

@@ -68,7 +68,7 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...;
         df = OnceDifferentiable(f!, vec(u0), vec(resid); autodiff)
     end
 
-    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
+    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,

ext/NonlinearSolveSIAMFANLEquationsExt.jl

@@ -2,51 +2,58 @@ module NonlinearSolveSIAMFANLEquationsExt
 
 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 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, arg
     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

src/NonlinearSolve.jl

@@ -237,7 +237,8 @@ export RadiusUpdateSchemes
 export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
     Broyden, Klement, LimitedMemoryBroyden
 export LeastSquaresOptimJL,
-    FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, FixedPointAccelerationJL, SpeedMappingJL, SIAMFANLEquationsJL
+    FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, FixedPointAccelerationJL, SpeedMappingJL,
+    SIAMFANLEquationsJL
 export NonlinearSolvePolyAlgorithm,
     RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg
 

src/extension_algs.jl

@@ -247,7 +247,6 @@ function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace =
 end
 
 """
-
     SpeedMappingJL(; σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false,
         orders::Vector{Int} = [3, 3, 2], time_limit::Real = 1000)
 
@@ -364,13 +363,11 @@ function FixedPointAccelerationJL(; algorithm = :Anderson, m = missing,
 end
 
 """
-
-    SIAMFANLEquationsJL(; method = :newton, autodiff = :central, show_trace = false, delta = 1e-3, linsolve = nothing)
+    SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
 
 ### Keyword Arguments
 
   - `method`: the choice of method for solving the nonlinear system.
-  - `show_trace`: whether to show the trace.
   - `delta`: initial pseudo time step, default is 1e-3.
   - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`.
 
@@ -380,16 +377,16 @@ end
   - `:pseudotransient`: Pseudo transient method.
   - `:secant`: Secant method for scalar equations.
 """
-@concrete struct SIAMFANLEquationsJL <: AbstractNonlinearAlgorithm
+@concrete struct SIAMFANLEquationsJL{L <: Union{Symbol, Nothing}} <:
+                 AbstractNonlinearSolveAlgorithm
     method::Symbol
-    show_trace::Bool
     delta
-    linsolve::Union{Symbol, Nothing}
+    linsolve::L
 end
 
-function SIAMFANLEquationsJL(; method = :newton, show_trace = false, delta = 1e-3, linsolve = nothing)
+function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
     if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing
         error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded")
     end
-  return SIAMFANLEquationsJL(method, show_trace, delta, linsolve)
+    return SIAMFANLEquationsJL(method, show_trace, delta, linsolve)
 end

test/siamfanlequations.jl

@@ -74,7 +74,7 @@ end
 f_tol(u, p) = u^2 - 2
 prob_tol = NonlinearProblem(f_tol, 1.0)
 for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-11]
-    for method = [:newton, :pseudotransient, :secant]
+    for method in [:newton, :pseudotransient, :secant]
         sol = solve(prob_tol, SIAMFANLEquationsJL(method = method), abstol = tol)
         @test abs(sol.u[1] - sqrt(2)) < tol
     end
@@ -141,12 +141,12 @@ f = NonlinearFunction(f!, jac = j!)
 p = A
 
 ProbN = NonlinearProblem(f, init, p)
-for method = [:newton, :pseudotransient]
+for method in [:newton, :pseudotransient]
     sol = solve(ProbN, SIAMFANLEquationsJL(method = method), reltol = 1e-8, abstol = 1e-8)
 end
 
 #= doesn't support complex numbers handling
 init = ones(Complex{Float64}, 152);
 ProbN = NonlinearProblem(f, init, p)
 sol = solve(ProbN, SIAMFANLEquationsJL(), reltol = 1e-8, abstol = 1e-8)
-=#
+=#