diff --git a/Project.toml b/Project.toml index f817aba9f..bd31721c1 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ BracketingNonlinearSolve = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e" CommonSolve = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2" ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471" DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" +DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63" FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" @@ -62,10 +63,12 @@ Aqua = "0.8" ArrayInterface = "7.16" BandedMatrices = "1.5" BenchmarkTools = "1.4" +BracketingNonlinearSolve = "1" CUDA = "5.5" CommonSolve = "0.2.4" ConcreteStructs = "0.2.3" DiffEqBase = "6.155.3" +DifferentiationInterface = "0.6.18" Enzyme = "0.13.11" ExplicitImports = "1.5" FastClosures = "0.3.2" @@ -87,6 +90,9 @@ NLsolve = "4.5" NaNMath = "1" NonlinearProblemLibrary = "0.1.2" NonlinearSolveBase = "1" +NonlinearSolveFirstOrder = "1" +NonlinearSolveQuasiNewton = "1" +NonlinearSolveSpectralMethods = "1" OrdinaryDiffEqTsit5 = "1.1.0" PETSc = "0.2" Pkg = "1.10" diff --git a/common/common_nlls_testing.jl b/common/common_nlls_testing.jl index 2888c56e5..86c33e8ac 100644 --- a/common/common_nlls_testing.jl +++ b/common/common_nlls_testing.jl @@ -47,3 +47,4 @@ prob_iip_vjp = NonlinearLeastSquaresProblem( ) export prob_oop, prob_iip, prob_oop_vjp, prob_iip_vjp +export true_function, θ_true, x, y_target, loss_function, θ_init diff --git a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl index e772043b3..a851c3a92 100644 --- a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl +++ b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl @@ -1,72 +1,84 @@ module NonlinearSolveFastLevenbergMarquardtExt -using ArrayInterface: ArrayInterface using FastClosures: @closure + +using ArrayInterface: ArrayInterface using FastLevenbergMarquardt: FastLevenbergMarquardt -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance -using NonlinearSolve: NonlinearSolve, FastLevenbergMarquardtJL -using SciMLBase: SciMLBase, NonlinearLeastSquaresProblem, NonlinearProblem, ReturnCode using StaticArraysCore: SArray -const FastLM = FastLevenbergMarquardt +using NonlinearSolveBase: NonlinearSolveBase +using NonlinearSolve: NonlinearSolve, FastLevenbergMarquardtJL +using SciMLBase: SciMLBase, AbstractNonlinearProblem, ReturnCode -@inline function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolve} - if linsolve === :cholesky - return FastLM.CholeskySolver(ArrayInterface.undefmatrix(x)) - elseif linsolve === :qr - return FastLM.QRSolver(eltype(x), length(x)) - else - throw(ArgumentError("Unknown FastLevenbergMarquardt Linear Solver: $linsolve")) - end -end -@inline _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, ::SArray) where {linsolve} = linsolve +const FastLM = FastLevenbergMarquardt -function SciMLBase.__solve(prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, - alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = nothing, - reltol = nothing, maxiters = 1000, termination_condition = nothing, kwargs...) - NonlinearSolve.__test_termination_condition( - termination_condition, :FastLevenbergMarquardt) +function SciMLBase.__solve( + prob::AbstractNonlinearProblem, alg::FastLevenbergMarquardtJL, args...; + alias_u0 = false, abstol = nothing, reltol = nothing, maxiters = 1000, + termination_condition = nothing, kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) - fn, u, resid = NonlinearSolve.__construct_extension_f( - prob; alias_u0, can_handle_oop = Val(prob.u0 isa SArray)) + f_wrapped, u, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0, can_handle_oop = Val(prob.u0 isa SArray) + ) f = if prob.u0 isa SArray - @closure (u, p) -> fn(u) + @closure (u, p) -> f_wrapped(u) else - @closure (du, u, p) -> fn(du, u) + @closure (du, u, p) -> f_wrapped(du, u) end - abstol = get_tolerance(abstol, eltype(u)) - reltol = get_tolerance(reltol, eltype(u)) - _jac_fn = NonlinearSolve.__construct_extension_jac( - prob, alg, u, resid; alg.autodiff, can_handle_oop = Val(prob.u0 isa SArray)) + abstol = NonlinearSolveBase.get_tolerance(abstol, eltype(u)) + reltol = NonlinearSolveBase.get_tolerance(reltol, eltype(u)) + + jac_fn_wrapped = NonlinearSolveBase.construct_extension_jac( + prob, alg, u, resid; alg.autodiff, can_handle_oop = Val(prob.u0 isa SArray) + ) jac_fn = if prob.u0 isa SArray - @closure (u, p) -> _jac_fn(u) + @closure (u, p) -> jac_fn_wrapped(u) else - @closure (J, u, p) -> _jac_fn(J, u) + @closure (J, u, p) -> jac_fn_wrapped(J, u) end - solver_kwargs = (; xtol = reltol, ftol = reltol, gtol = abstol, maxit = maxiters, + solver_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) + alg.maxscale, alg.factorupdate, alg.minfactor, alg.maxfactor + ) if prob.u0 isa SArray res, fx, info, iter, nfev, njev = FastLM.lmsolve( - f, jac_fn, prob.u0; solver_kwargs...) + f, jac_fn, prob.u0; solver_kwargs... + ) LM, solver = nothing, nothing else J = prob.f.jac_prototype === nothing ? similar(u, length(resid), length(u)) : zero(prob.f.jac_prototype) - solver = _fast_lm_solver(alg, u) + + solver = if alg.linsolve === :cholesky + FastLM.CholeskySolver(ArrayInterface.undefmatrix(u)) + elseif alg.linsolve === :qr + FastLM.QRSolver(eltype(u), length(u)) + else + throw(ArgumentError("Unknown FastLevenbergMarquardt Linear Solver: \ + $(Meta.quot(alg.linsolve))")) + end + LM = FastLM.LMWorkspace(u, resid, J) res, fx, info, iter, nfev, njev, LM, solver = FastLM.lmsolve!( - f, jac_fn, LM; solver, solver_kwargs...) + f, jac_fn, LM; solver, solver_kwargs... + ) end stats = SciMLBase.NLStats(nfev, njev, -1, -1, iter) retcode = info == -1 ? ReturnCode.MaxIters : ReturnCode.Success - return SciMLBase.build_solution(prob, alg, res, fx; retcode, - original = (res, fx, info, iter, nfev, njev, LM, solver), stats) + return SciMLBase.build_solution( + prob, alg, res, fx; retcode, + original = (res, fx, info, iter, nfev, njev, LM, solver), stats + ) end end diff --git a/ext/NonlinearSolveFixedPointAccelerationExt.jl b/ext/NonlinearSolveFixedPointAccelerationExt.jl index 2d36d7f56..6ff9f240c 100644 --- a/ext/NonlinearSolveFixedPointAccelerationExt.jl +++ b/ext/NonlinearSolveFixedPointAccelerationExt.jl @@ -1,41 +1,51 @@ module NonlinearSolveFixedPointAccelerationExt -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance +using FixedPointAcceleration: FixedPointAcceleration, fixed_point + +using NonlinearSolveBase: NonlinearSolveBase using NonlinearSolve: NonlinearSolve, FixedPointAccelerationJL using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode -using FixedPointAcceleration: FixedPointAcceleration, fixed_point -function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL, args...; +function SciMLBase.__solve( + prob::NonlinearProblem, alg::FixedPointAccelerationJL, args...; abstol = nothing, maxiters = 1000, alias_u0::Bool = false, - show_trace::Val{PrintReports} = Val(false), - termination_condition = nothing, kwargs...) where {PrintReports} - NonlinearSolve.__test_termination_condition( - termination_condition, :FixedPointAccelerationJL) + show_trace::Val = Val(false), termination_condition = nothing, kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) + + f, u0, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0, make_fixed_point = Val(true), force_oop = Val(true) + ) - f, u0, resid = NonlinearSolve.__construct_extension_f( - prob; alias_u0, make_fixed_point = Val(true), force_oop = Val(true)) - tol = get_tolerance(abstol, eltype(u0)) + tol = NonlinearSolveBase.get_tolerance(abstol, eltype(u0)) - sol = fixed_point(f, u0; Algorithm = alg.algorithm, MaxIter = maxiters, MaxM = alg.m, + sol = fixed_point( + f, u0; Algorithm = alg.algorithm, MaxIter = maxiters, MaxM = alg.m, ConvergenceMetricThreshold = tol, ExtrapolationPeriod = alg.extrapolation_period, - Dampening = alg.dampening, PrintReports, ReplaceInvalids = alg.replace_invalids, - ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true) + Dampening = alg.dampening, PrintReports = show_trace isa Val{true}, + ReplaceInvalids = alg.replace_invalids, + ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true + ) if sol.FixedPoint_ === missing u0 = prob.u0 isa Number ? u0[1] : u0 - resid = NonlinearSolve.evaluate_f(prob, u0) + resid = NonlinearSolveBase.Utils.evaluate_f(prob, u0) res = u0 converged = false else res = prob.u0 isa Number ? first(sol.FixedPoint_) : reshape(sol.FixedPoint_, size(prob.u0)) - resid = NonlinearSolve.evaluate_f(prob, res) + resid = NonlinearSolveBase.Utils.evaluate_f(prob, res) converged = maximum(abs, resid) ≤ tol end - return SciMLBase.build_solution(prob, alg, res, resid; original = sol, + return SciMLBase.build_solution( + prob, alg, res, resid; original = sol, retcode = converged ? ReturnCode.Success : ReturnCode.Failure, - stats = SciMLBase.NLStats(sol.Iterations_, 0, 0, 0, sol.Iterations_)) + stats = SciMLBase.NLStats(sol.Iterations_, 0, 0, 0, sol.Iterations_) + ) end end diff --git a/ext/NonlinearSolveLeastSquaresOptimExt.jl b/ext/NonlinearSolveLeastSquaresOptimExt.jl index 51563eb83..0df265963 100644 --- a/ext/NonlinearSolveLeastSquaresOptimExt.jl +++ b/ext/NonlinearSolveLeastSquaresOptimExt.jl @@ -1,79 +1,70 @@ module NonlinearSolveLeastSquaresOptimExt -using ConcreteStructs: @concrete using LeastSquaresOptim: LeastSquaresOptim -using NonlinearSolveBase: NonlinearSolveBase, TraceMinimal, get_tolerance + +using NonlinearSolveBase: NonlinearSolveBase, TraceMinimal using NonlinearSolve: NonlinearSolve, LeastSquaresOptimJL -using SciMLBase: SciMLBase, NonlinearLeastSquaresProblem, NonlinearProblem, ReturnCode +using SciMLBase: SciMLBase, AbstractNonlinearProblem, ReturnCode const LSO = LeastSquaresOptim -@inline function _lso_solver(::LeastSquaresOptimJL{alg, ls}) where {alg, ls} - linsolve = ls === :qr ? LSO.QR() : - (ls === :cholesky ? LSO.Cholesky() : (ls === :lsmr ? LSO.LSMR() : nothing)) - if alg === :lm - return LSO.LevenbergMarquardt(linsolve) - elseif alg === :dogleg - return LSO.Dogleg(linsolve) - else - throw(ArgumentError("Unknown LeastSquaresOptim Algorithm: $alg")) - end -end - -@concrete struct LeastSquaresOptimJLCache - prob - alg - allocated_prob - kwargs -end - -function Base.show(io::IO, cache::LeastSquaresOptimJLCache) - print(io, "LeastSquaresOptimJLCache()") -end - -function SciMLBase.reinit!(cache::LeastSquaresOptimJLCache, args...; kwargs...) - error("Reinitialization not supported for LeastSquaresOptimJL.") -end - -function SciMLBase.__init(prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, - alg::LeastSquaresOptimJL, args...; alias_u0 = false, abstol = nothing, - show_trace::Val{ShT} = Val(false), trace_level = TraceMinimal(), - reltol = nothing, store_trace::Val{StT} = Val(false), maxiters = 1000, - termination_condition = nothing, kwargs...) where {ShT, StT} - NonlinearSolve.__test_termination_condition(termination_condition, :LeastSquaresOptim) +function SciMLBase.__solve( + prob::AbstractNonlinearProblem, alg::LeastSquaresOptimJL, args...; + alias_u0 = false, abstol = nothing, reltol = nothing, maxiters = 1000, + trace_level = TraceMinimal(), termination_condition = nothing, + show_trace::Val = Val(false), store_trace::Val = Val(false), kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) - f!, u, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) - abstol = get_tolerance(abstol, eltype(u)) - reltol = get_tolerance(reltol, eltype(u)) + f!, u, resid = NonlinearSolveBase.construct_extension_function_wrapper(prob; alias_u0) + abstol = NonlinearSolveBase.get_tolerance(abstol, eltype(u)) + reltol = NonlinearSolveBase.get_tolerance(reltol, eltype(u)) if prob.f.jac === nothing && alg.autodiff isa Symbol - lsoprob = LSO.LeastSquaresProblem(; x = u, f!, y = resid, alg.autodiff, - J = prob.f.jac_prototype, output_length = length(resid)) + lsoprob = LSO.LeastSquaresProblem(; + x = u, f!, y = resid, alg.autodiff, J = prob.f.jac_prototype, + output_length = length(resid) + ) else - g! = NonlinearSolve.__construct_extension_jac(prob, alg, u, resid; alg.autodiff) + g! = NonlinearSolveBase.construct_extension_jac(prob, alg, u, resid; alg.autodiff) lsoprob = LSO.LeastSquaresProblem(; x = u, f!, y = resid, g!, J = prob.f.jac_prototype, - output_length = length(resid)) + output_length = length(resid) + ) end - allocated_prob = LSO.LeastSquaresProblemAllocated(lsoprob, _lso_solver(alg)) + linsolve = alg.ls === :qr ? LSO.QR() : + (alg.ls === :cholesky ? LSO.Cholesky() : + (alg.ls === :lsmr ? LSO.LSMR() : nothing)) - return LeastSquaresOptimJLCache(prob, - alg, - allocated_prob, - (; x_tol = reltol, f_tol = abstol, g_tol = abstol, iterations = maxiters, - show_trace = ShT, store_trace = StT, show_every = trace_level.print_frequency)) -end + lso_solver = if alg.alg === :lm + LSO.LevenbergMarquardt(linsolve) + elseif alg.alg === :dogleg + LSO.Dogleg(linsolve) + else + throw(ArgumentError("Unknown LeastSquaresOptim Algorithm: $(Meta.quot(alg.alg))")) + end + + allocated_prob = LSO.LeastSquaresProblemAllocated(lsoprob, lso_solver(alg)) + res = LSO.optimize!( + allocated_prob; + x_tol = reltol, f_tol = abstol, g_tol = abstol, iterations = maxiters, + show_trace = show_trace isa Val{true}, store_trace = store_trace isa Val{true}, + show_every = trace_level.print_frequency + ) -function SciMLBase.solve!(cache::LeastSquaresOptimJLCache) - res = LSO.optimize!(cache.allocated_prob; cache.kwargs...) - maxiters = cache.kwargs[:iterations] retcode = res.x_converged || res.f_converged || res.g_converged ? ReturnCode.Success : (res.iterations ≥ maxiters ? ReturnCode.MaxIters : ReturnCode.ConvergenceFailure) stats = SciMLBase.NLStats(res.f_calls, res.g_calls, -1, -1, res.iterations) + + f!(resid, res.minimizer) + return SciMLBase.build_solution( - cache.prob, cache.alg, res.minimizer, res.ssr / 2; retcode, original = res, stats) + prob, alg, res.minimizer, resid; retcode, original = res, stats + ) end end diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl index 88adf5753..fc0f919a1 100644 --- a/ext/NonlinearSolveMINPACKExt.jl +++ b/ext/NonlinearSolveMINPACKExt.jl @@ -1,52 +1,67 @@ module NonlinearSolveMINPACKExt +using FastClosures: @closure using MINPACK: MINPACK -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance + +using NonlinearSolveBase: NonlinearSolveBase using NonlinearSolve: NonlinearSolve, CMINPACK using SciMLBase: SciMLBase, NonlinearLeastSquaresProblem, NonlinearProblem, ReturnCode -using FastClosures: @closure function SciMLBase.__solve( - prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, alg::CMINPACK, - args...; 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 {ShT, StT} - NonlinearSolve.__test_termination_condition(termination_condition, :CMINPACK) - - _f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) - f! = @closure (du, u) -> (_f!(du, u); Cint(0)) + prob::Union{NonlinearLeastSquaresProblem, NonlinearProblem}, alg::CMINPACK, args...; + abstol = nothing, maxiters = 1000, alias_u0::Bool = false, + show_trace::Val = Val(false), store_trace::Val = Val(false), + termination_condition = nothing, kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) + + f_wrapped!, u0, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0 + ) + resid_size = size(resid) + f! = @closure (du, u) -> (f_wrapped!(du, u); Cint(0)) m = length(resid) - method = ifelse(alg.method === :auto, - ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hybr), alg.method) + method = ifelse( + alg.method === :auto, + ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hybr), alg.method + ) - show_trace = ShT - tracing = StT - tol = get_tolerance(abstol, eltype(u0)) + show_trace = show_trace isa Val{true} + tracing = store_trace isa Val{true} + tol = NonlinearSolveBase.get_tolerance(abstol, eltype(u0)) if alg.autodiff === missing && prob.f.jac === nothing original = MINPACK.fsolve( - f!, u0, m; tol, show_trace, tracing, method, iterations = maxiters) + f!, u0, m; tol, show_trace, tracing, method, iterations = maxiters + ) else autodiff = alg.autodiff === missing ? nothing : alg.autodiff - _jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; autodiff) - jac! = @closure (J, u) -> (_jac!(J, u); Cint(0)) + jac_wrapped! = NonlinearSolveBase.construct_extension_jac( + prob, alg, u0, resid; autodiff + ) + jac! = @closure (J, u) -> (jac_wrapped!(J, u); Cint(0)) original = MINPACK.fsolve( - f!, jac!, u0, m; tol, show_trace, tracing, method, iterations = maxiters) + f!, jac!, u0, m; tol, show_trace, tracing, method, iterations = maxiters + ) end u = original.x - resid_ = original.f - objective = maximum(abs, resid_) + resid = original.f + objective = maximum(abs, resid) retcode = ifelse(objective ≤ tol, ReturnCode.Success, ReturnCode.Failure) # These are only meaningful if `store_trace = Val(true)` - stats = SciMLBase.NLStats(original.trace.f_calls, original.trace.g_calls, - original.trace.g_calls, original.trace.g_calls, -1) + stats = SciMLBase.NLStats( + original.trace.f_calls, original.trace.g_calls, + original.trace.g_calls, original.trace.g_calls, -1 + ) - u_ = prob.u0 isa Number ? original.x[1] : reshape(original.x, size(prob.u0)) - resid_ = prob.u0 isa Number ? resid_[1] : reshape(resid_, size(resid)) - return SciMLBase.build_solution(prob, alg, u_, resid_; retcode, original, stats) + u = prob.u0 isa Number ? original.x[1] : reshape(original.x, size(prob.u0)) + resid = prob.u0 isa Number ? resid[1] : reshape(resid, resid_size) + return SciMLBase.build_solution(prob, alg, u, resid; retcode, original, stats) end end diff --git a/ext/NonlinearSolveNLSolversExt.jl b/ext/NonlinearSolveNLSolversExt.jl index a9f41c87c..95aa0c8c3 100644 --- a/ext/NonlinearSolveNLSolversExt.jl +++ b/ext/NonlinearSolveNLSolversExt.jl @@ -1,74 +1,63 @@ module NonlinearSolveNLSolversExt using ADTypes: ADTypes, AutoFiniteDiff, AutoForwardDiff, AutoPolyesterForwardDiff +using DifferentiationInterface: DifferentiationInterface, Constant using FastClosures: @closure -using FiniteDiff: FiniteDiff -using ForwardDiff: ForwardDiff -using LinearAlgebra: norm + using NLSolvers: NLSolvers, NEqOptions, NEqProblem -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance + +using NonlinearSolveBase: NonlinearSolveBase using NonlinearSolve: NonlinearSolve, NLSolversJL using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode -function SciMLBase.__solve(prob::NonlinearProblem, alg::NLSolversJL, args...; - abstol = nothing, reltol = nothing, maxiters = 1000, - alias_u0::Bool = false, termination_condition = nothing, kwargs...) - NonlinearSolve.__test_termination_condition(termination_condition, :NLSolversJL) +const DI = DifferentiationInterface + +function SciMLBase.__solve( + prob::NonlinearProblem, alg::NLSolversJL, args...; + abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0::Bool = false, + termination_condition = nothing, kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) - abstol = get_tolerance(abstol, eltype(prob.u0)) - reltol = get_tolerance(reltol, eltype(prob.u0)) + abstol = NonlinearSolveBase.get_tolerance(abstol, eltype(prob.u0)) + reltol = NonlinearSolveBase.get_tolerance(reltol, eltype(prob.u0)) options = NEqOptions(; maxiter = maxiters, f_abstol = abstol, f_reltol = reltol) if prob.u0 isa Number - f_scalar = @closure x -> prob.f(x, prob.p) - - if alg.autodiff === nothing - if ForwardDiff.can_dual(typeof(prob.u0)) - autodiff_concrete = :forwarddiff - else - autodiff_concrete = :finitediff - end - else - if alg.autodiff isa AutoForwardDiff || alg.autodiff isa AutoPolyesterForwardDiff - autodiff_concrete = :forwarddiff - elseif alg.autodiff isa AutoFiniteDiff - autodiff_concrete = :finitediff - else - error("Only ForwardDiff or FiniteDiff autodiff is supported.") - end - end + f_scalar = Base.Fix2(prob.f, prob.p) + autodiff = NonlinearSolveBase.select_jacobian_autodiff(prob, alg.autodiff) + + prep = DifferentiationInterface.prepare_derivative( + prob.f, autodiff, prob.u0, Constant(prob.p) + ) - if autodiff_concrete === :forwarddiff - fj_scalar = @closure (Jx, x) -> begin - T = typeof(ForwardDiff.Tag(prob.f, eltype(x))) - x_dual = ForwardDiff.Dual{T}(x, one(x)) - y = prob.f(x_dual, prob.p) - return ForwardDiff.value(y), ForwardDiff.extract_derivative(T, y) - end - else - fj_scalar = @closure (Jx, x) -> begin - _f = Base.Fix2(prob.f, prob.p) - return _f(x), FiniteDiff.finite_difference_derivative(_f, x) - end + fj_scalar = @closure (Jx, x) -> begin + return DifferentiationInterface.value_and_derivative( + prob.f, prep, autodiff, x, Constant(prob.p) + ) end prob_obj = NLSolvers.ScalarObjective(; f = f_scalar, fg = fj_scalar) prob_nlsolver = NEqProblem(prob_obj; inplace = false) res = NLSolvers.solve(prob_nlsolver, prob.u0, alg.method, options) - retcode = ifelse(norm(res.info.best_residual, Inf) ≤ abstol, - ReturnCode.Success, ReturnCode.MaxIters) + retcode = ifelse( + maximum(abs, res.info.best_residual) ≤ abstol, + ReturnCode.Success, ReturnCode.MaxIters + ) stats = SciMLBase.NLStats(-1, -1, -1, -1, res.info.iter) return SciMLBase.build_solution( prob, alg, res.info.solution, res.info.best_residual; - retcode, original = res, stats) + retcode, original = res, stats + ) end - f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) - - jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; alg.autodiff) + f!, u0, resid = NonlinearSolveBase.construct_extension_function_wrapper(prob; alias_u0) + jac! = NonlinearSolveBase.construct_extension_jac(prob, alg, u0, resid; alg.autodiff) FJ_vector! = @closure (Fx, Jx, x) -> begin f!(Fx, x) @@ -82,11 +71,15 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLSolversJL, args...; res = NLSolvers.solve(prob_nlsolver, u0, alg.method, options) retcode = ifelse( - norm(res.info.best_residual, Inf) ≤ abstol, ReturnCode.Success, ReturnCode.MaxIters) + maximum(abs, res.info.best_residual) ≤ abstol, + ReturnCode.Success, ReturnCode.MaxIters + ) stats = SciMLBase.NLStats(-1, -1, -1, -1, res.info.iter) - return SciMLBase.build_solution(prob, alg, res.info.solution, res.info.best_residual; - retcode, original = res, stats) + return SciMLBase.build_solution( + prob, alg, res.info.solution, res.info.best_residual; + retcode, original = res, stats + ) end end diff --git a/ext/NonlinearSolveNLsolveExt.jl b/ext/NonlinearSolveNLsolveExt.jl index 49c1f337a..1b0beb993 100644 --- a/ext/NonlinearSolveNLsolveExt.jl +++ b/ext/NonlinearSolveNLsolveExt.jl @@ -1,44 +1,51 @@ module NonlinearSolveNLsolveExt using LineSearches: Static -using NonlinearSolveBase: NonlinearSolveBase, TraceMinimal, get_tolerance -using NonlinearSolve: NonlinearSolve, NLsolveJL using NLsolve: NLsolve, OnceDifferentiable, nlsolve + +using NonlinearSolveBase: NonlinearSolveBase, Utils, TraceMinimal +using NonlinearSolve: NonlinearSolve, NLsolveJL using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode function SciMLBase.__solve( - prob::NonlinearProblem, alg::NLsolveJL, args...; abstol = nothing, - maxiters = 1000, alias_u0::Bool = false, termination_condition = nothing, - store_trace::Val{StT} = Val(false), show_trace::Val{ShT} = Val(false), - trace_level = TraceMinimal(), kwargs...) where {StT, ShT} - NonlinearSolve.__test_termination_condition(termination_condition, :NLsolveJL) + prob::NonlinearProblem, alg::NLsolveJL, args...; + abstol = nothing, maxiters = 1000, alias_u0::Bool = false, + termination_condition = nothing, trace_level = TraceMinimal(), + store_trace::Val = Val(false), show_trace::Val = Val(false), kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) - f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) + f!, u0, resid = NonlinearSolveBase.construct_extension_function_wrapper(prob; alias_u0) if prob.f.jac === nothing && alg.autodiff isa Symbol df = OnceDifferentiable(f!, u0, resid; alg.autodiff) else autodiff = alg.autodiff isa Symbol ? nothing : alg.autodiff - jac! = NonlinearSolve.__construct_extension_jac(prob, alg, u0, resid; autodiff) + jac! = NonlinearSolveBase.construct_extension_jac(prob, alg, u0, resid; autodiff) if prob.f.jac_prototype === nothing J = similar( - u0, promote_type(eltype(u0), eltype(resid)), length(u0), length(resid)) + u0, promote_type(eltype(u0), eltype(resid)), length(u0), length(resid) + ) else J = zero(prob.f.jac_prototype) end - df = OnceDifferentiable(f!, jac!, vec(u0), vec(resid), J) + df = OnceDifferentiable(f!, jac!, Utils.safe_vec(u0), Utils.safe_vec(resid), J) end - abstol = get_tolerance(abstol, eltype(u0)) - show_trace = ShT - store_trace = StT + abstol = NonlinearSolveBase.get_tolerance(abstol, eltype(u0)) + show_trace = show_trace isa Val{true} + store_trace = store_trace isa Val{true} extended_trace = !(trace_level.trace_mode isa Val{:minimal}) linesearch = alg.linesearch === missing ? Static() : alg.linesearch - original = nlsolve(df, vec(u0); ftol = abstol, iterations = maxiters, alg.method, - store_trace, extended_trace, linesearch, alg.linsolve, alg.factor, - alg.autoscale, alg.m, alg.beta, show_trace) + original = nlsolve( + df, vec(u0); + ftol = abstol, iterations = maxiters, alg.method, store_trace, extended_trace, + linesearch, alg.linsolve, alg.factor, alg.autoscale, alg.m, alg.beta, show_trace + ) f!(vec(resid), original.zero) u = prob.u0 isa Number ? original.zero[1] : reshape(original.zero, size(prob.u0)) @@ -46,8 +53,10 @@ function SciMLBase.__solve( retcode = original.x_converged || original.f_converged ? ReturnCode.Success : ReturnCode.Failure - stats = SciMLBase.NLStats(original.f_calls, original.g_calls, original.g_calls, - original.g_calls, original.iterations) + stats = SciMLBase.NLStats( + original.f_calls, original.g_calls, original.g_calls, + original.g_calls, original.iterations + ) return SciMLBase.build_solution(prob, alg, u, resid; retcode, original, stats) end diff --git a/ext/NonlinearSolvePETScExt.jl b/ext/NonlinearSolvePETScExt.jl index 9146c0b61..7b36cecce 100644 --- a/ext/NonlinearSolvePETScExt.jl +++ b/ext/NonlinearSolvePETScExt.jl @@ -1,23 +1,31 @@ module NonlinearSolvePETScExt using FastClosures: @closure + using MPI: MPI -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance -using NonlinearSolve: NonlinearSolve, PETScSNES using PETSc: PETSc + +using NonlinearSolveBase: NonlinearSolveBase +using NonlinearSolve: NonlinearSolve, PETScSNES using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode + using SparseArrays: AbstractSparseMatrix function SciMLBase.__solve( - prob::NonlinearProblem, alg::PETScSNES, args...; abstol = nothing, reltol = nothing, + prob::NonlinearProblem, alg::PETScSNES, args...; + abstol = nothing, reltol = nothing, maxiters = 1000, alias_u0::Bool = false, termination_condition = nothing, - show_trace::Val{ShT} = Val(false), kwargs...) where {ShT} + show_trace::Val = Val(false), kwargs... +) # XXX: https://petsc.org/release/manualpages/SNES/SNESSetConvergenceTest/ - termination_condition === nothing || - error("`PETScSNES` does not support termination conditions!") - - _f!, u0, resid = NonlinearSolve.__construct_extension_f(prob; alias_u0) - T = eltype(prob.u0) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg; abs_norm_supported = false + ) + + f_wrapped!, u0, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0 + ) + T = eltype(u0) @assert T ∈ PETSc.scalar_types if alg.petsclib === missing @@ -35,8 +43,8 @@ function SciMLBase.__solve( end PETSc.initialized(petsclib) || PETSc.initialize(petsclib) - abstol = get_tolerance(abstol, T) - reltol = get_tolerance(reltol, T) + abstol = NonlinearSolveBase.get_tolerance(abstol, T) + reltol = NonlinearSolveBase.get_tolerance(reltol, T) nf = Ref{Int}(0) @@ -44,77 +52,82 @@ function SciMLBase.__solve( nf[] += 1 fx = cfx isa Ptr{Nothing} ? PETSc.unsafe_localarray(T, cfx; read = false) : cfx x = cx isa Ptr{Nothing} ? PETSc.unsafe_localarray(T, cx; write = false) : cx - _f!(fx, x) + f_wrapped!(fx, x) Base.finalize(fx) Base.finalize(x) return end - snes = PETSc.SNES{T}(petsclib, + snes = PETSc.SNES{T}( + petsclib, alg.mpi_comm === missing ? MPI.COMM_SELF : alg.mpi_comm; - alg.snes_options..., snes_monitor = ShT, snes_rtol = reltol, - snes_atol = abstol, snes_max_it = maxiters) + alg.snes_options..., snes_monitor = show_trace isa Val{true}, snes_rtol = reltol, + snes_atol = abstol, snes_max_it = maxiters + ) PETSc.setfunction!(snes, f!, PETSc.VecSeq(zero(u0))) - if alg.autodiff === missing && prob.f.jac === nothing - _jac! = nothing - njac = Ref{Int}(-1) - else + njac = Ref{Int}(-1) + if alg.autodiff !== missing || prob.f.jac !== nothing autodiff = alg.autodiff === missing ? nothing : alg.autodiff if prob.u0 isa Number - _jac! = NonlinearSolve.__construct_extension_jac( - prob, alg, prob.u0, prob.u0; autodiff) + jac! = NonlinearSolveBase.construct_extension_jac( + prob, alg, prob.u0, prob.u0; autodiff + ) J_init = zeros(T, 1, 1) else - _jac!, J_init = NonlinearSolve.__construct_extension_jac( - prob, alg, u0, resid; autodiff, initial_jacobian = Val(true)) + jac!, J_init = NonlinearSolveBase.construct_extension_jac( + prob, alg, u0, resid; autodiff, initial_jacobian = Val(true) + ) end njac = Ref{Int}(0) if J_init isa AbstractSparseMatrix PJ = PETSc.MatSeqAIJ(J_init) - jac! = @closure (cx, J, _, user_ctx) -> begin + jac_fn! = @closure (cx, J, _, user_ctx) -> begin njac[] += 1 x = cx isa Ptr{Nothing} ? PETSc.unsafe_localarray(T, cx; write = false) : cx if J isa PETSc.AbstractMat - _jac!(user_ctx.jacobian, x) + jac!(user_ctx.jacobian, x) copyto!(J, user_ctx.jacobian) PETSc.assemble(J) else - _jac!(J, x) + jac!(J, x) end Base.finalize(x) return end - PETSc.setjacobian!(snes, jac!, PJ, PJ) + PETSc.setjacobian!(snes, jac_fn!, PJ, PJ) snes.user_ctx = (; jacobian = J_init) else PJ = PETSc.MatSeqDense(J_init) - jac! = @closure (cx, J, _, user_ctx) -> begin + jac_fn! = @closure (cx, J, _, user_ctx) -> begin njac[] += 1 x = cx isa Ptr{Nothing} ? PETSc.unsafe_localarray(T, cx; write = false) : cx - _jac!(J, x) + jac!(J, x) Base.finalize(x) J isa PETSc.AbstractMat && PETSc.assemble(J) return end - PETSc.setjacobian!(snes, jac!, PJ, PJ) + PETSc.setjacobian!(snes, jac_fn!, PJ, PJ) end end res = PETSc.solve!(u0, snes) - _f!(resid, res) - u_ = prob.u0 isa Number ? res[1] : res - resid_ = prob.u0 isa Number ? resid[1] : resid + f_wrapped!(resid, res) + u_res = prob.u0 isa Number ? res[1] : res + resid_res = prob.u0 isa Number ? resid[1] : resid objective = maximum(abs, resid) # XXX: Return Code from PETSc retcode = ifelse(objective ≤ abstol, ReturnCode.Success, ReturnCode.Failure) - return SciMLBase.build_solution(prob, alg, u_, resid_; retcode, original = snes, - stats = SciMLBase.NLStats(nf[], njac[], -1, -1, -1)) + return SciMLBase.build_solution( + prob, alg, u_res, resid_res; + retcode, original = snes, + stats = SciMLBase.NLStats(nf[], njac[], -1, -1, -1) + ) end end diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl index 2468064fb..bfd2bd72d 100644 --- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl +++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl @@ -1,13 +1,14 @@ module NonlinearSolveSIAMFANLEquationsExt using FastClosures: @closure -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance -using NonlinearSolve: NonlinearSolve, SIAMFANLEquationsJL -using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode using SIAMFANLEquations: SIAMFANLEquations, aasol, nsol, nsoli, nsolsc, ptcsol, ptcsoli, ptcsolsc, secant -@inline function __siam_fanl_equations_retcode_mapping(sol) +using NonlinearSolveBase: NonlinearSolveBase +using NonlinearSolve: NonlinearSolve, SIAMFANLEquationsJL +using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode + +function siamfanlequations_retcode_mapping(sol) if sol.errcode == 0 return ReturnCode.Success elseif sol.errcode == 10 @@ -16,102 +17,129 @@ using SIAMFANLEquations: SIAMFANLEquations, aasol, nsol, nsoli, nsolsc, ptcsol, return ReturnCode.Failure elseif sol.errcode == -1 return ReturnCode.Default + else + error("Unknown SIAMFANLEquations return code: $(sol.errcode)") end end -@inline function __zeros_like(x, args...) +function zeros_like(x, args...) z = similar(x, args...) - fill!(z, zero(eltype(x))) + fill!(z, false) return z 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) +function siamfanlequations_stats_mapping(method, sol) ((method === :pseudotransient) || (method === :anderson)) && return nothing return SciMLBase.NLStats( - sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)) + sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm) + ) end -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} - NonlinearSolve.__test_termination_condition(termination_condition, :SIAMFANLEquationsJL) +function SciMLBase.__solve( + prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...; + abstol = nothing, reltol = nothing, alias_u0::Bool = false, maxiters = 1000, + termination_condition = nothing, show_trace = Val(false), kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) (; method, delta, linsolve, m, beta) = alg T = eltype(prob.u0) - atol = get_tolerance(abstol, T) - rtol = get_tolerance(reltol, T) + atol = NonlinearSolveBase.get_tolerance(abstol, T) + rtol = NonlinearSolveBase.get_tolerance(reltol, T) + + printerr = show_trace isa Val{true} if prob.u0 isa Number - f = @closure u -> prob.f(u, prob.p) + f = Base.Fix2(prob.f, prob.p) if method == :newton - sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT) + sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr) elseif method == :pseudotransient sol = ptcsolsc( - f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol, printerr = ShT) + f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol, printerr + ) elseif method == :secant - sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT) + sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr) elseif method == :anderson - f_aa, u, _ = NonlinearSolve.__construct_extension_f( - prob; alias_u0, make_fixed_point = Val(true)) - sol = aasol(f_aa, u, m, __zeros_like(u, 1, 2 * m + 4); - maxit = maxiters, atol, rtol, beta) + f_aa, u, _ = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0, make_fixed_point = Val(true) + ) + sol = aasol( + f_aa, u, m, zeros_like(u, 1, 2 * m + 4); + maxit = maxiters, atol, rtol, beta + ) end else - f, u, resid = NonlinearSolve.__construct_extension_f( - prob; alias_u0, make_fixed_point = Val(method == :anderson)) + f, u, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0, make_fixed_point = Val(method == :anderson) + ) N = length(u) - FS = __zeros_like(u, N) + FS = zeros_like(u, N) # Jacobian Free Newton Krylov if linsolve !== nothing # Allocate ahead for Krylov basis - JVS = linsolve == :gmres ? __zeros_like(u, N, 3) : __zeros_like(u, N) + JVS = linsolve == :gmres ? zeros_like(u, N, 3) : zeros_like(u, N) linsolve_alg = String(linsolve) if method == :newton - sol = nsoli(f, u, FS, JVS; lsolver = linsolve_alg, - maxit = maxiters, atol, rtol, printerr = ShT) + sol = nsoli( + f, u, FS, JVS; lsolver = linsolve_alg, + maxit = maxiters, atol, rtol, printerr + ) elseif method == :pseudotransient - sol = ptcsoli(f, u, FS, JVS; lsolver = linsolve_alg, - maxit = maxiters, atol, rtol, printerr = ShT) + sol = ptcsoli( + f, u, FS, JVS; lsolver = linsolve_alg, + maxit = maxiters, atol, rtol, printerr + ) end else if prob.f.jac === nothing && alg.autodiff === missing - FPS = __zeros_like(u, N, N) + FPS = zeros_like(u, N, N) if method == :newton - sol = nsol(f, u, FS, FPS; sham = 1, atol, rtol, - maxit = maxiters, printerr = ShT) + sol = nsol( + f, u, FS, FPS; sham = 1, atol, rtol, maxit = maxiters, printerr + ) elseif method == :pseudotransient - sol = ptcsol(f, u, FS, FPS; atol, rtol, maxit = maxiters, - delta0 = delta, printerr = ShT) + sol = ptcsol( + f, u, FS, FPS; + atol, rtol, maxit = maxiters, delta0 = delta, printerr + ) elseif method == :anderson sol = aasol( - f, u, m, zeros(T, N, 2 * m + 4); atol, rtol, maxit = maxiters, beta) + f, u, m, zeros(T, N, 2 * m + 4); + atol, rtol, maxit = maxiters, beta + ) end else autodiff = alg.autodiff === missing ? nothing : alg.autodiff FPS = prob.f.jac_prototype !== nothing ? zero(prob.f.jac_prototype) : - __zeros_like(u, N, N) - jac = NonlinearSolve.__construct_extension_jac( - prob, alg, u, resid; autodiff) + zeros_like(u, N, N) + jac = NonlinearSolveBase.construct_extension_jac( + prob, alg, u, resid; autodiff + ) AJ! = @closure (J, u, x) -> jac(J, x) if method == :newton - sol = nsol(f, u, FS, FPS, AJ!; sham = 1, atol, - rtol, maxit = maxiters, printerr = ShT) + sol = nsol( + f, u, FS, FPS, AJ!; sham = 1, atol, + rtol, maxit = maxiters, printerr + ) elseif method == :pseudotransient - sol = ptcsol(f, u, FS, FPS, AJ!; atol, rtol, maxit = maxiters, - delta0 = delta, printerr = ShT) + sol = ptcsol( + f, u, FS, FPS, AJ!; atol, rtol, maxit = maxiters, + delta0 = delta, printerr + ) end end end end - retcode = __siam_fanl_equations_retcode_mapping(sol) - stats = __siam_fanl_equations_stats_mapping(method, sol) + retcode = siamfanlequations_retcode_mapping(sol) + stats = siamfanlequations_stats_mapping(method, sol) res = prob.u0 isa Number && method === :anderson ? sol.solution[1] : sol.solution - resid = NonlinearSolve.evaluate_f(prob, res) + resid = NonlinearSolveBase.Utils.evaluate_f(prob, res) return SciMLBase.build_solution(prob, alg, res, resid; retcode, stats, original = sol) end diff --git a/ext/NonlinearSolveSpeedMappingExt.jl b/ext/NonlinearSolveSpeedMappingExt.jl index ff9b4683b..c0f39607b 100644 --- a/ext/NonlinearSolveSpeedMappingExt.jl +++ b/ext/NonlinearSolveSpeedMappingExt.jl @@ -1,30 +1,41 @@ module NonlinearSolveSpeedMappingExt -using NonlinearSolveBase: NonlinearSolveBase, get_tolerance +using SpeedMapping: speedmapping + +using NonlinearSolveBase: NonlinearSolveBase using NonlinearSolve: NonlinearSolve, SpeedMappingJL using SciMLBase: SciMLBase, NonlinearProblem, ReturnCode -using SpeedMapping: speedmapping -function SciMLBase.__solve(prob::NonlinearProblem, alg::SpeedMappingJL, args...; +function SciMLBase.__solve( + prob::NonlinearProblem, alg::SpeedMappingJL, args...; abstol = nothing, maxiters = 1000, alias_u0::Bool = false, - maxtime = nothing, store_trace::Val{store_info} = Val(false), - termination_condition = nothing, kwargs...) where {store_info} - NonlinearSolve.__test_termination_condition(termination_condition, :SpeedMappingJL) + maxtime = nothing, store_trace::Val = Val(false), + termination_condition = nothing, kwargs... +) + NonlinearSolveBase.assert_extension_supported_termination_condition( + termination_condition, alg + ) - m!, u, resid = NonlinearSolve.__construct_extension_f( - prob; alias_u0, make_fixed_point = Val(true)) - tol = get_tolerance(abstol, eltype(u)) + m!, u, resid = NonlinearSolveBase.construct_extension_function_wrapper( + prob; alias_u0, make_fixed_point = Val(true) + ) + tol = NonlinearSolveBase.get_tolerance(abstol, eltype(u)) time_limit = ifelse(maxtime === nothing, 1000, maxtime) - sol = speedmapping(u; m!, tol, Lp = Inf, maps_limit = maxiters, alg.orders, - alg.check_obj, store_info, alg.σ_min, alg.stabilize, time_limit) + sol = speedmapping( + u; m!, tol, Lp = Inf, maps_limit = maxiters, alg.orders, + alg.check_obj, store_info = store_trace isa Val{true}, alg.σ_min, alg.stabilize, + time_limit + ) res = prob.u0 isa Number ? first(sol.minimizer) : sol.minimizer - resid = NonlinearSolve.evaluate_f(prob, res) + resid = NonlinearSolveBase.Utils.evaluate_f(prob, res) - return SciMLBase.build_solution(prob, alg, res, resid; original = sol, - retcode = sol.converged ? ReturnCode.Success : ReturnCode.Failure, - stats = SciMLBase.NLStats(sol.maps, 0, 0, 0, sol.maps)) + return SciMLBase.build_solution( + prob, alg, res, resid; + original = sol, stats = SciMLBase.NLStats(sol.maps, 0, 0, 0, sol.maps), + retcode = ifelse(sol.converged, ReturnCode.Success, ReturnCode.Failure) + ) end end diff --git a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl index a3ac4614d..a08384677 100644 --- a/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl +++ b/lib/NonlinearSolveBase/src/NonlinearSolveBase.jl @@ -46,6 +46,7 @@ include("jacobian.jl") include("linear_solve.jl") include("timer_outputs.jl") include("tracing.jl") +include("wrappers.jl") include("descent/common.jl") include("descent/newton.jl") @@ -66,6 +67,9 @@ include("solve.jl") @compat(public, (InternalAPI, supports_line_search, supports_trust_region, set_du!)) @compat(public, (construct_linear_solver, needs_square_A, needs_concrete_A)) @compat(public, (construct_jacobian_cache,)) +@compat(public, + (assert_extension_supported_termination_condition, + construct_extension_function_wrapper, construct_extension_jac)) export TraceMinimal, TraceWithJacobianConditionNumber, TraceAll diff --git a/lib/NonlinearSolveBase/src/wrappers.jl b/lib/NonlinearSolveBase/src/wrappers.jl new file mode 100644 index 000000000..0b96836ad --- /dev/null +++ b/lib/NonlinearSolveBase/src/wrappers.jl @@ -0,0 +1,113 @@ +function assert_extension_supported_termination_condition( + termination_condition, alg; abs_norm_supported = true +) + no_termination_condition = termination_condition === nothing + no_termination_condition && return nothing + abs_norm_supported && termination_condition isa AbsNormTerminationMode && return nothing + throw(AssertionError("`$(nameof(typeof(alg)))` does not support termination conditions!")) +end + +function construct_extension_function_wrapper( + prob::AbstractNonlinearProblem; alias_u0::Bool = false, + can_handle_oop::Val = Val(false), can_handle_scalar::Val = Val(false), + make_fixed_point::Val = Val(false), force_oop::Val = Val(false) +) + if can_handle_oop isa Val{false} && can_handle_scalar isa Val{true} + error("Incorrect Specification: OOP not supported but scalar supported.") + end + + resid = Utils.evaluate_f(prob, prob.u0) + u0 = can_handle_scalar isa Val{true} || !(prob.u0 isa Number) ? + Utils.maybe_unaliased(prob.u0, alias_u0) : [prob.u0] + + fₚ = if make_fixed_point isa Val{true} + if SciMLBase.isinplace(prob) + @closure (du, u) -> begin + prob.f(du, u, prob.p) + du .+= u + return du + end + else + @closure u -> prob.f(u, prob.p) .+ u + end + else + if SciMLBase.isinplace(prob) + @closure (du, u) -> begin + prob.f(du, u, prob.p) + return du + end + else + Base.Fix2(prob.f, prob.p) + end + end + + f_flat_structure = if SciMLBase.isinplace(prob) + u0_size, du_size = size(u0), size(resid) + @closure (du, u) -> begin + fₚ(reshape(du, du_size), reshape(u, u0_size)) + return du + end + else + if prob.u0 isa Number + if can_handle_scalar isa Val{true} + fₚ + elseif can_handle_oop isa Val{true} + @closure u -> [fₚ(first(u))] + else + @closure (du, u) -> begin + du[1] = fₚ(first(u)) + return du + end + end + else + u0_size = size(u0) + if can_handle_oop isa Val{true} + @closure u -> vec(fₚ(reshape(u, u0_size))) + else + @closure (du, u) -> begin + copyto!(du, fₚ(reshape(u, u0_size))) + return du + end + end + end + end + + f_final = if force_oop isa Val{true} && applicable(f_flat_structure, u0, u0) + resid = resid isa Number ? [resid] : Utils.safe_vec(resid) + du = Utils.safe_vec(zero(resid)) + @closure u -> begin + f_flat_structure(du, u) + return du + end + else + f_flat_structure + end + + return f_final, Utils.safe_vec(u0), (resid isa Number ? [resid] : Utils.safe_vec(resid)) +end + +function construct_extension_jac( + prob, alg, u0, fu; + can_handle_oop::Val = Val(false), can_handle_scalar::Val = Val(false), + autodiff = nothing, initial_jacobian = Val(false), kwargs... +) + autodiff = select_jacobian_autodiff(prob, autodiff) + + Jₚ = construct_jacobian_cache( + prob, alg, prob.f, fu, u0, prob.p; + stats = NLStats(0, 0, 0, 0, 0), autodiff, kwargs... + ) + + J_no_scalar = can_handle_scalar isa Val{false} && prob.u0 isa Number ? + @closure(u->[Jₚ(u[1])]) : Jₚ + + J_final = if can_handle_oop isa Val{false} && !SciMLBase.isinplace(prob) + @closure((J, u)->copyto!(J, J_no_scalar(u))) + else + J_no_scalar + end + + initial_jacobian isa Val{false} && return J_final + + return J_final, Jₚ(nothing) +end diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index 87bc2a5b3..cc71ef03d 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -49,7 +49,7 @@ const SII = SymbolicIndexingInterface include("helpers.jl") include("polyalg.jl") -# include("extension_algs.jl") +include("extension_algs.jl") include("default.jl") diff --git a/src/extension_algs.jl b/src/extension_algs.jl index fb922fb14..1fad79839 100644 --- a/src/extension_algs.jl +++ b/src/extension_algs.jl @@ -1,469 +1,490 @@ -# # This file only include the algorithm struct to be exported by NonlinearSolve.jl. The main -# # functionality is implemented as package extensions -# """ -# LeastSquaresOptimJL(alg = :lm; linsolve = nothing, autodiff::Symbol = :central) - -# Wrapper over [LeastSquaresOptim.jl](https://github.com/matthieugomez/LeastSquaresOptim.jl) -# for solving `NonlinearLeastSquaresProblem`. - -# ### Arguments - -# - `alg`: Algorithm to use. Can be `:lm` or `:dogleg`. - -# ### Keyword Arguments - -# - `linsolve`: Linear solver to use. Can be `:qr`, `:cholesky` or `:lsmr`. If `nothing`, -# then `LeastSquaresOptim.jl` will choose the best linear solver based on the Jacobian -# structure. -# - `autodiff`: Automatic differentiation / Finite Differences. Can be `:central` or -# `:forward`. - -# !!! note - -# This algorithm is only available if `LeastSquaresOptim.jl` is installed and loaded. -# """ -# struct LeastSquaresOptimJL{alg, linsolve} <: AbstractNonlinearSolveExtensionAlgorithm -# autodiff -# end - -# function LeastSquaresOptimJL(alg = :lm; linsolve = nothing, autodiff = :central) -# @assert alg in (:lm, :dogleg) -# @assert linsolve === nothing || linsolve in (:qr, :cholesky, :lsmr) -# autodiff isa Symbol && @assert autodiff in (:central, :forward) - -# if Base.get_extension(@__MODULE__, :NonlinearSolveLeastSquaresOptimExt) === nothing -# error("LeastSquaresOptimJL requires LeastSquaresOptim.jl to be loaded") -# end - -# return LeastSquaresOptimJL{alg, linsolve}(autodiff) -# end - -# """ -# FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6, -# factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt, -# minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, -# autodiff = nothing) - -# Wrapper over [FastLevenbergMarquardt.jl](https://github.com/kamesy/FastLevenbergMarquardt.jl) -# for solving `NonlinearLeastSquaresProblem`. For details about the other keyword arguments -# see the documentation for `FastLevenbergMarquardt.jl`. - -# !!! warning - -# This is not really the fastest solver. It is called that since the original package -# is called "Fast". `LevenbergMarquardt()` is almost always a better choice. - -# ### Arguments - -# - `linsolve`: Linear solver to use. Can be `:qr` or `:cholesky`. - -# ### Keyword Arguments - -# - `autodiff`: determines the backend used for the Jacobian. Note that this argument is -# ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to -# `nothing` which means that a default is selected according to the problem specification! - -# !!! note - -# This algorithm is only available if `FastLevenbergMarquardt.jl` is installed and loaded. -# """ -# @concrete struct FastLevenbergMarquardtJL{linsolve} <: -# AbstractNonlinearSolveExtensionAlgorithm -# autodiff -# factor -# factoraccept -# factorreject -# factorupdate::Symbol -# minscale -# maxscale -# minfactor -# maxfactor -# end - -# function FastLevenbergMarquardtJL( -# linsolve::Symbol = :cholesky; factor = 1e-6, factoraccept = 13.0, -# factorreject = 3.0, factorupdate = :marquardt, minscale = 1e-12, -# maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, autodiff = nothing) -# @assert linsolve in (:qr, :cholesky) -# @assert factorupdate in (:marquardt, :nielson) - -# if Base.get_extension(@__MODULE__, :NonlinearSolveFastLevenbergMarquardtExt) === nothing -# error("FastLevenbergMarquardtJL requires FastLevenbergMarquardt.jl to be loaded") -# end - -# return FastLevenbergMarquardtJL{linsolve}(autodiff, factor, factoraccept, factorreject, -# factorupdate, minscale, maxscale, minfactor, maxfactor) -# end - -# """ -# CMINPACK(; method::Symbol = :auto, autodiff = missing) - -# ### Keyword Arguments - -# - `method`: the choice of method for the solver. -# - `autodiff`: Defaults to `missing`, which means we will default to letting `MINPACK` -# construct the jacobian if `f.jac` is not provided. In other cases, we use it to generate -# a jacobian similar to other NonlinearSolve solvers. - -# ### Submethod Choice - -# The keyword argument `method` can take on different value depending on which method of -# `fsolve` you are calling. The standard choices of `method` are: - -# - `:hybr`: Modified version of Powell's algorithm. Uses MINPACK routine -# [`hybrd1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd1.c) -# - `:lm`: Levenberg-Marquardt. Uses MINPACK routine -# [`lmdif1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif1.c) -# - `:lmdif`: Advanced Levenberg-Marquardt (more options available with `; kwargs...`). See -# MINPACK routine [`lmdif`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif.c) -# for more information -# - `:hybrd`: Advanced modified version of Powell's algorithm (more options available with -# `; kwargs...`). See MINPACK routine -# [`hybrd`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd.c) -# for more information - -# If a Jacobian is supplied as part of the [`NonlinearFunction`](@ref nonlinearfunctions), -# then the following methods are allowed: - -# - `:hybr`: Advanced modified version of Powell's algorithm with user supplied Jacobian. -# Additional arguments are available via `; kwargs...`. See MINPACK routine -# [`hybrj`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrj.c) -# for more information -# - `:lm`: Advanced Levenberg-Marquardt with user supplied Jacobian. Additional arguments -# are available via `; kwargs...`. See MINPACK routine -# [`lmder`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmder.c) -# for more information - -# The default choice of `:auto` selects `:hybr` for NonlinearProblem and `:lm` for -# NonlinearLeastSquaresProblem. - -# !!! note - -# This algorithm is only available if `MINPACK.jl` is installed and loaded. -# """ -# @concrete struct CMINPACK <: AbstractNonlinearSolveExtensionAlgorithm -# method::Symbol -# autodiff -# end - -# function CMINPACK(; method::Symbol = :auto, autodiff = missing) -# if Base.get_extension(@__MODULE__, :NonlinearSolveMINPACKExt) === nothing -# error("CMINPACK requires MINPACK.jl to be loaded") -# end -# return CMINPACK(method, autodiff) -# end - -# """ -# NLsolveJL(; method = :trust_region, autodiff = :central, linesearch = Static(), -# linsolve = (x, A, b) -> copyto!(x, A\\b), factor = one(Float64), autoscale = true, -# m = 10, beta = one(Float64)) - -# ### Keyword Arguments - -# - `method`: the choice of method for solving the nonlinear system. -# - `autodiff`: the choice of method for generating the Jacobian. Defaults to `:central` or -# central differencing via FiniteDiff.jl. The other choices are `:forward` or `ADTypes` -# similar to other solvers in NonlinearSolve. -# - `linesearch`: the line search method to be used within the solver method. The choices -# are line search types from -# [LineSearches.jl](https://github.com/JuliaNLSolvers/LineSearches.jl). -# - `linsolve`: a function `linsolve(x, A, b)` that solves `Ax = b`. -# - `factor`: determines the size of the initial trust region. This size is set to the -# product of factor and the euclidean norm of `u0` if nonzero, or else to factor itself. -# - `autoscale`: if true, then the variables will be automatically rescaled. The scaling -# factors are the norms of the Jacobian columns. -# - `m`: the amount of history in the Anderson method. Naive "Picard"-style iteration can be -# achieved by setting m=0, but that isn't advisable for contractions whose Lipschitz -# constants are close to 1. If convergence fails, though, you may consider lowering it. -# - `beta`: It is also known as DIIS or Pulay mixing, this method is based on the -# acceleration of the fixed-point iteration xₙ₊₁ = xₙ + beta*f(xₙ), where by default -# beta = 1. - -# !!! warning - -# Line Search Algorithms from [`LineSearch.jl`](https://github.com/SciML/LineSearch.jl) -# aren't supported by `NLsolveJL`. Instead, use the line search algorithms from -# [`LineSearches.jl`](https://github.com/JuliaNLSolvers/LineSearches.jl). - -# ### Submethod Choice - -# Choices for methods in `NLsolveJL`: - -# - `:anderson`: Anderson-accelerated fixed-point iteration -# - `:broyden`: Broyden's quasi-Newton method -# - `:newton`: Classical Newton method with an optional line search -# - `:trust_region`: Trust region Newton method (the default choice) - -# For more information on these arguments, consult the -# [NLsolve.jl documentation](https://github.com/JuliaNLSolvers/NLsolve.jl). - -# !!! note - -# This algorithm is only available if `NLsolve.jl` is installed and loaded. -# """ -# @concrete struct NLsolveJL <: AbstractNonlinearSolveExtensionAlgorithm -# method::Symbol -# autodiff -# linesearch -# linsolve -# factor -# autoscale::Bool -# m::Int -# beta -# end - -# function NLsolveJL(; method = :trust_region, autodiff = :central, linesearch = missing, -# linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, -# autoscale = true, m = 10, beta = one(Float64)) -# if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolveExt) === nothing -# error("NLsolveJL requires NLsolve.jl to be loaded") -# end - -# if autodiff isa Symbol && autodiff !== :central && autodiff !== :forward -# error("`autodiff` must be `:central` or `:forward`.") -# end - -# return NLsolveJL(method, autodiff, linesearch, linsolve, factor, autoscale, m, beta) -# end - -# """ -# NLSolversJL(method; autodiff = nothing) -# NLSolversJL(; method, autodiff = nothing) - -# Wrapper over NLSolvers.jl Nonlinear Equation Solvers. We automatically construct the -# jacobian function and supply it to the solver. - -# ### Arguments - -# - `method`: the choice of method for solving the nonlinear system. See the documentation -# for NLSolvers.jl for more information. -# - `autodiff`: the choice of method for generating the Jacobian. Defaults to `nothing` -# which means that a default is selected according to the problem specification. Can be -# any valid ADTypes.jl autodiff type (conditional on that backend being supported in -# NonlinearSolve.jl). -# """ -# struct NLSolversJL{M, AD} <: AbstractNonlinearSolveExtensionAlgorithm -# method::M -# autodiff::AD - -# function NLSolversJL(method, autodiff) -# if Base.get_extension(@__MODULE__, :NonlinearSolveNLSolversExt) === nothing -# error("NLSolversJL requires NLSolvers.jl to be loaded") -# end - -# return new{typeof(method), typeof(autodiff)}(method, autodiff) -# end -# end - -# NLSolversJL(method; autodiff = nothing) = NLSolversJL(method, autodiff) -# NLSolversJL(; method, autodiff = nothing) = NLSolversJL(method, autodiff) - -# """ -# SpeedMappingJL(; σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false, -# orders::Vector{Int} = [3, 3, 2], time_limit::Real = 1000) - -# Wrapper over [SpeedMapping.jl](https://nicolasl-s.github.io/SpeedMapping.jl/) for solving -# Fixed Point Problems. We allow using this algorithm to solve root finding problems as well. - -# ### Keyword Arguments - -# - `σ_min`: Setting to `1` may avoid stalling (see [lepage2021alternating](@cite)). -# - `stabilize`: performs a stabilization mapping before extrapolating. Setting to `true` -# may improve the performance for applications like accelerating the EM or MM algorithms -# (see [lepage2021alternating](@cite)). -# - `check_obj`: In case of NaN or Inf values, the algorithm restarts at the best past -# iterate. -# - `orders`: determines ACX's alternating order. Must be between `1` and `3` (where `1` -# means no extrapolation). The two recommended orders are `[3, 2]` and `[3, 3, 2]`, the -# latter being potentially better for highly non-linear applications (see -# [lepage2021alternating](@cite)). -# - `time_limit`: time limit for the algorithm. - -# !!! note - -# This algorithm is only available if `SpeedMapping.jl` is installed and loaded. -# """ -# @concrete struct SpeedMappingJL <: AbstractNonlinearSolveExtensionAlgorithm -# σ_min -# stabilize::Bool -# check_obj::Bool -# orders::Vector{Int} -# end - -# function SpeedMappingJL(; σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false, -# orders::Vector{Int} = [3, 3, 2]) -# if Base.get_extension(@__MODULE__, :NonlinearSolveSpeedMappingExt) === nothing -# error("SpeedMappingJL requires SpeedMapping.jl to be loaded") -# end - -# return SpeedMappingJL(σ_min, stabilize, check_obj, orders) -# end - -# """ -# FixedPointAccelerationJL(; algorithm = :Anderson, m = missing, -# condition_number_threshold = missing, extrapolation_period = missing, -# replace_invalids = :NoAction) - -# Wrapper over [FixedPointAcceleration.jl](https://s-baumann.github.io/FixedPointAcceleration.jl/) -# for solving Fixed Point Problems. We allow using this algorithm to solve root finding -# problems as well. - -# ### Keyword Arguments - -# - `algorithm`: The algorithm to use. Can be `:Anderson`, `:MPE`, `:RRE`, `:VEA`, `:SEA`, -# `:Simple`, `:Aitken` or `:Newton`. -# - `m`: The number of previous iterates to use for the extrapolation. Only valid for -# `:Anderson`. -# - `condition_number_threshold`: The condition number threshold for Least Squares Problem. -# Only valid for `:Anderson`. -# - `extrapolation_period`: The number of iterates between extrapolations. Only valid for -# `:MPE`, `:RRE`, `:VEA` and `:SEA`. Defaults to `7` for `:MPE` & `:RRE`, and `6` for -# `:SEA` and `:VEA`. For `:SEA` and `:VEA`, this must be a multiple of `2`. -# - `replace_invalids`: The method to use for replacing invalid iterates. Can be -# `:ReplaceInvalids`, `:ReplaceVector` or `:NoAction`. - -# !!! note - -# This algorithm is only available if `FixedPointAcceleration.jl` is installed and loaded. -# """ -# @concrete struct FixedPointAccelerationJL <: AbstractNonlinearSolveExtensionAlgorithm -# algorithm::Symbol -# extrapolation_period::Int -# replace_invalids::Symbol -# dampening -# m::Int -# condition_number_threshold -# end - -# function FixedPointAccelerationJL(; -# algorithm = :Anderson, m = missing, condition_number_threshold = missing, -# extrapolation_period = missing, replace_invalids = :NoAction, dampening = 1.0) -# if Base.get_extension(@__MODULE__, :NonlinearSolveFixedPointAccelerationExt) === nothing -# error("FixedPointAccelerationJL requires FixedPointAcceleration.jl to be loaded") -# end - -# @assert algorithm in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) -# @assert replace_invalids in (:ReplaceInvalids, :ReplaceVector, :NoAction) - -# if algorithm !== :Anderson -# if condition_number_threshold !== missing -# error("`condition_number_threshold` is only valid for Anderson acceleration") -# end -# if m !== missing -# error("`m` is only valid for Anderson acceleration") -# end -# end -# condition_number_threshold === missing && (condition_number_threshold = 1e3) -# m === missing && (m = 10) - -# if algorithm !== :MPE && algorithm !== :RRE && algorithm !== :VEA && algorithm !== :SEA -# if extrapolation_period !== missing -# error("`extrapolation_period` is only valid for MPE, RRE, VEA and SEA") -# end -# end -# if extrapolation_period === missing -# extrapolation_period = algorithm === :SEA || algorithm === :VEA ? 6 : 7 -# else -# if (algorithm === :SEA || algorithm === :VEA) && extrapolation_period % 2 != 0 -# error("`extrapolation_period` must be multiples of 2 for SEA and VEA") -# end -# end - -# return FixedPointAccelerationJL(algorithm, extrapolation_period, replace_invalids, -# dampening, m, condition_number_threshold) -# end - -# """ -# SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, -# autodiff = missing) - -# ### Keyword Arguments - -# - `method`: the choice of method for solving the nonlinear system. -# - `delta`: initial pseudo time step, default is 1e-3. -# - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`. -# - `m`: Depth for Anderson acceleration, default as 0 for Picard iteration. -# - `beta`: Anderson mixing parameter, change f(x) to (1-beta)x+beta*f(x), -# equivalent to accelerating damped Picard iteration. -# - `autodiff`: Defaults to `missing`, which means we will default to letting -# `SIAMFANLEquations` construct the jacobian if `f.jac` is not provided. In other cases, -# we use it to generate a jacobian similar to other NonlinearSolve solvers. - -# ### Submethod Choice - -# - `:newton`: Classical Newton method. -# - `:pseudotransient`: Pseudo transient method. -# - `:secant`: Secant method for scalar equations. -# - `:anderson`: Anderson acceleration for fixed point iterations. - -# !!! note - -# This algorithm is only available if `SIAMFANLEquations.jl` is installed and loaded. -# """ -# @concrete struct SIAMFANLEquationsJL <: AbstractNonlinearSolveExtensionAlgorithm -# method::Symbol -# delta -# linsolve <: Union{Symbol, Nothing} -# m::Int -# beta -# autodiff -# end - -# function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing, -# m = 0, beta = 1.0, autodiff = missing) -# if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing -# error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded") -# end -# return SIAMFANLEquationsJL(method, delta, linsolve, m, beta, autodiff) -# end - -# """ -# PETScSNES(; petsclib = missing, autodiff = nothing, mpi_comm = missing, kwargs...) - -# Wrapper over [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl) SNES solvers. - -# ### Keyword Arguments - -# - `petsclib`: PETSc library to use. If set to `missing`, then we will use the first -# available PETSc library in `PETSc.petsclibs` based on the problem element type. -# - `autodiff`: the choice of method for generating the Jacobian. Defaults to `nothing` -# which means that a default is selected according to the problem specification. Can be -# any valid ADTypes.jl autodiff type (conditional on that backend being supported in -# NonlinearSolve.jl). If set to `missing`, then PETSc computes the Jacobian using finite -# differences. -# - `mpi_comm`: MPI communicator to use. If set to `missing`, then we will use -# `MPI.COMM_SELF`. -# - `kwargs`: Keyword arguments to be passed to the PETSc SNES solver. See [PETSc SNES -# documentation](https://petsc.org/release/manual/snes/) and -# [SNESSetFromOptions](https://petsc.org/release/manualpages/SNES/SNESSetFromOptions) -# for more information. - -# ### Options via `CommonSolve.solve` - -# These options are forwarded from `solve` to the PETSc SNES solver. If these are provided to -# `kwargs`, then they will be ignored. - -# | `solve` option | PETSc SNES option | -# |:-------------- |:----------------- | -# | `atol` | `snes_atol` | -# | `rtol` | `snes_rtol` | -# | `maxiters` | `snes_max_it` | -# | `show_trace` | `snes_monitor` | - -# !!! note - -# This algorithm is only available if `PETSc.jl` is installed and loaded. -# """ -# @concrete struct PETScSNES <: AbstractNonlinearSolveExtensionAlgorithm -# petsclib -# mpi_comm -# autodiff <: Union{Missing, Nothing, ADTypes.AbstractADType} -# snes_options -# end - -# function PETScSNES(; petsclib = missing, autodiff = nothing, mpi_comm = missing, kwargs...) -# if Base.get_extension(@__MODULE__, :NonlinearSolvePETScExt) === nothing -# error("PETScSNES requires PETSc.jl to be loaded") -# end -# return PETScSNES(petsclib, mpi_comm, autodiff, kwargs) -# end +# This file only include the algorithm struct to be exported by NonlinearSolve.jl. The main +# functionality is implemented as package extensions +""" + LeastSquaresOptimJL(alg = :lm; linsolve = nothing, autodiff::Symbol = :central) + +Wrapper over [LeastSquaresOptim.jl](https://github.com/matthieugomez/LeastSquaresOptim.jl) +for solving `NonlinearLeastSquaresProblem`. + +### Arguments + + - `alg`: Algorithm to use. Can be `:lm` or `:dogleg`. + +### Keyword Arguments + + - `linsolve`: Linear solver to use. Can be `:qr`, `:cholesky` or `:lsmr`. If `nothing`, + then `LeastSquaresOptim.jl` will choose the best linear solver based on the Jacobian + structure. + - `autodiff`: Automatic differentiation / Finite Differences. Can be `:central` or + `:forward`. + +!!! note + + This algorithm is only available if `LeastSquaresOptim.jl` is installed and loaded. +""" +@concrete struct LeastSquaresOptimJL <: AbstractNonlinearSolveAlgorithm + autodiff + alg::Symbol + linsolve <: Union{Symbol, Nothing} + name::Symbol +end + +function LeastSquaresOptimJL(alg = :lm; linsolve = nothing, autodiff = :central) + @assert alg in (:lm, :dogleg) + @assert linsolve === nothing || linsolve in (:qr, :cholesky, :lsmr) + autodiff isa Symbol && @assert autodiff in (:central, :forward) + + if Base.get_extension(@__MODULE__, :NonlinearSolveLeastSquaresOptimExt) === nothing + error("`LeastSquaresOptimJL` requires `LeastSquaresOptim.jl` to be loaded") + end + + return LeastSquaresOptimJL(autodiff, alg, linsolve, :LeastSquaresOptimJL) +end + +""" + FastLevenbergMarquardtJL( + linsolve::Symbol = :cholesky; + factor = 1e-6, factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt, + minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, + autodiff = nothing + ) + +Wrapper over [FastLevenbergMarquardt.jl](https://github.com/kamesy/FastLevenbergMarquardt.jl) +for solving `NonlinearLeastSquaresProblem`. For details about the other keyword arguments +see the documentation for `FastLevenbergMarquardt.jl`. + +### Arguments + + - `linsolve`: Linear solver to use. Can be `:qr` or `:cholesky`. + +### Keyword Arguments + + - `autodiff`: determines the backend used for the Jacobian. Note that this argument is + ignored if an analytical Jacobian is passed, as that will be used instead. Defaults to + `nothing` which means that a default is selected according to the problem specification! + +!!! note + + This algorithm is only available if `FastLevenbergMarquardt.jl` is installed and loaded. +""" +@concrete struct FastLevenbergMarquardtJL <: AbstractNonlinearSolveAlgorithm + autodiff + linsolve::Symbol + factor + factoraccept + factorreject + factorupdate::Symbol + minscale + maxscale + minfactor + maxfactor + name::Symbol +end + +function FastLevenbergMarquardtJL( + linsolve::Symbol = :cholesky; factor = 1e-6, factoraccept = 13.0, + factorreject = 3.0, factorupdate = :marquardt, minscale = 1e-12, + maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32, autodiff = nothing +) + @assert linsolve in (:qr, :cholesky) + @assert factorupdate in (:marquardt, :nielson) + + if Base.get_extension(@__MODULE__, :NonlinearSolveFastLevenbergMarquardtExt) === nothing + error("`FastLevenbergMarquardtJL` requires `FastLevenbergMarquardt.jl` to be loaded") + end + + return FastLevenbergMarquardtJL( + autodiff, linsolve, factor, factoraccept, factorreject, + factorupdate, minscale, maxscale, minfactor, maxfactor, :FastLevenbergMarquardtJL + ) +end + +""" + CMINPACK(; method::Symbol = :auto, autodiff = missing) + +### Keyword Arguments + + - `method`: the choice of method for the solver. + - `autodiff`: Defaults to `missing`, which means we will default to letting `MINPACK` + construct the jacobian if `f.jac` is not provided. In other cases, we use it to generate + a jacobian similar to other NonlinearSolve solvers. + +### Submethod Choice + +The keyword argument `method` can take on different value depending on which method of +`fsolve` you are calling. The standard choices of `method` are: + + - `:hybr`: Modified version of Powell's algorithm. Uses MINPACK routine + [`hybrd1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd1.c) + - `:lm`: Levenberg-Marquardt. Uses MINPACK routine + [`lmdif1`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif1.c) + - `:lmdif`: Advanced Levenberg-Marquardt (more options available with `; kwargs...`). See + MINPACK routine [`lmdif`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmdif.c) + for more information + - `:hybrd`: Advanced modified version of Powell's algorithm (more options available with + `; kwargs...`). See MINPACK routine + [`hybrd`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrd.c) + for more information + +If a Jacobian is supplied as part of the [`NonlinearFunction`](@ref nonlinearfunctions), +then the following methods are allowed: + + - `:hybr`: Advanced modified version of Powell's algorithm with user supplied Jacobian. + Additional arguments are available via `; kwargs...`. See MINPACK routine + [`hybrj`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/hybrj.c) + for more information + - `:lm`: Advanced Levenberg-Marquardt with user supplied Jacobian. Additional arguments + are available via `; kwargs...`. See MINPACK routine + [`lmder`](https://github.com/devernay/cminpack/blob/d1f5f5a273862ca1bbcf58394e4ac060d9e22c76/lmder.c) + for more information + +The default choice of `:auto` selects `:hybr` for NonlinearProblem and `:lm` for +NonlinearLeastSquaresProblem. + +!!! note + + This algorithm is only available if `MINPACK.jl` is installed and loaded. +""" +@concrete struct CMINPACK <: AbstractNonlinearSolveAlgorithm + method::Symbol + autodiff + name::Symbol +end + +function CMINPACK(; method::Symbol = :auto, autodiff = missing) + if Base.get_extension(@__MODULE__, :NonlinearSolveMINPACKExt) === nothing + error("`CMINPACK` requires `MINPACK.jl` to be loaded") + end + return CMINPACK(method, autodiff, :CMINPACK) +end + +""" + NLsolveJL(; + method = :trust_region, autodiff = :central, linesearch = Static(), + linsolve = (x, A, b) -> copyto!(x, A\\b), factor = one(Float64), autoscale = true, + m = 10, beta = one(Float64) + ) + +### Keyword Arguments + + - `method`: the choice of method for solving the nonlinear system. + - `autodiff`: the choice of method for generating the Jacobian. Defaults to `:central` or + central differencing via FiniteDiff.jl. The other choices are `:forward` or `ADTypes` + similar to other solvers in NonlinearSolve. + - `linesearch`: the line search method to be used within the solver method. The choices + are line search types from + [LineSearches.jl](https://github.com/JuliaNLSolvers/LineSearches.jl). + - `linsolve`: a function `linsolve(x, A, b)` that solves `Ax = b`. + - `factor`: determines the size of the initial trust region. This size is set to the + product of factor and the euclidean norm of `u0` if nonzero, or else to factor itself. + - `autoscale`: if true, then the variables will be automatically rescaled. The scaling + factors are the norms of the Jacobian columns. + - `m`: the amount of history in the Anderson method. Naive "Picard"-style iteration can be + achieved by setting m=0, but that isn't advisable for contractions whose Lipschitz + constants are close to 1. If convergence fails, though, you may consider lowering it. + - `beta`: It is also known as DIIS or Pulay mixing, this method is based on the + acceleration of the fixed-point iteration xₙ₊₁ = xₙ + beta*f(xₙ), where by default + beta = 1. + +!!! warning + + Line Search Algorithms from [`LineSearch.jl`](https://github.com/SciML/LineSearch.jl) + aren't supported by `NLsolveJL`. Instead, use the line search algorithms from + [`LineSearches.jl`](https://github.com/JuliaNLSolvers/LineSearches.jl). + +### Submethod Choice + +Choices for methods in `NLsolveJL`: + + - `:anderson`: Anderson-accelerated fixed-point iteration + - `:broyden`: Broyden's quasi-Newton method + - `:newton`: Classical Newton method with an optional line search + - `:trust_region`: Trust region Newton method (the default choice) + +For more information on these arguments, consult the +[NLsolve.jl documentation](https://github.com/JuliaNLSolvers/NLsolve.jl). + +!!! note + + This algorithm is only available if `NLsolve.jl` is installed and loaded. +""" +@concrete struct NLsolveJL <: AbstractNonlinearSolveAlgorithm + method::Symbol + autodiff + linesearch + linsolve + factor + autoscale::Bool + m::Int + beta + name::Symbol +end + +function NLsolveJL(; + method = :trust_region, autodiff = :central, linesearch = missing, beta = 1.0, + linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, autoscale = true, m = 10 +) + if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolveExt) === nothing + error("`NLsolveJL` requires `NLsolve.jl` to be loaded") + end + + if autodiff isa Symbol && autodiff !== :central && autodiff !== :forward + error("`autodiff` must be `:central` or `:forward`.") + end + + return NLsolveJL( + method, autodiff, linesearch, linsolve, factor, autoscale, m, beta, :NLsolveJL + ) +end + +""" + NLSolversJL(method; autodiff = nothing) + NLSolversJL(; method, autodiff = nothing) + +Wrapper over NLSolvers.jl Nonlinear Equation Solvers. We automatically construct the +jacobian function and supply it to the solver. + +### Arguments + + - `method`: the choice of method for solving the nonlinear system. See the documentation + for NLSolvers.jl for more information. + - `autodiff`: the choice of method for generating the Jacobian. Defaults to `nothing` + which means that a default is selected according to the problem specification. Can be + any valid ADTypes.jl autodiff type (conditional on that backend being supported in + NonlinearSolve.jl). +""" +struct NLSolversJL{M, AD} <: AbstractNonlinearSolveAlgorithm + method::M + autodiff::AD + name::Symbol + + function NLSolversJL(method, autodiff) + if Base.get_extension(@__MODULE__, :NonlinearSolveNLSolversExt) === nothing + error("NLSolversJL requires NLSolvers.jl to be loaded") + end + + return new{typeof(method), typeof(autodiff)}(method, autodiff, :NLSolversJL) + end +end + +NLSolversJL(method; autodiff = nothing) = NLSolversJL(method, autodiff) +NLSolversJL(; method, autodiff = nothing) = NLSolversJL(method, autodiff) + +""" + SpeedMappingJL(; + σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false, + orders::Vector{Int} = [3, 3, 2] + ) + +Wrapper over [SpeedMapping.jl](https://nicolasl-s.github.io/SpeedMapping.jl/) for solving +Fixed Point Problems. We allow using this algorithm to solve root finding problems as well. + +### Keyword Arguments + + - `σ_min`: Setting to `1` may avoid stalling (see [lepage2021alternating](@cite)). + - `stabilize`: performs a stabilization mapping before extrapolating. Setting to `true` + may improve the performance for applications like accelerating the EM or MM algorithms + (see [lepage2021alternating](@cite)). + - `check_obj`: In case of NaN or Inf values, the algorithm restarts at the best past + iterate. + - `orders`: determines ACX's alternating order. Must be between `1` and `3` (where `1` + means no extrapolation). The two recommended orders are `[3, 2]` and `[3, 3, 2]`, the + latter being potentially better for highly non-linear applications (see + [lepage2021alternating](@cite)). + +!!! note + + This algorithm is only available if `SpeedMapping.jl` is installed and loaded. +""" +@concrete struct SpeedMappingJL <: AbstractNonlinearSolveAlgorithm + σ_min + stabilize::Bool + check_obj::Bool + orders::Vector{Int} + name::Symbol +end + +function SpeedMappingJL(; + σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false, + orders::Vector{Int} = [3, 3, 2] +) + if Base.get_extension(@__MODULE__, :NonlinearSolveSpeedMappingExt) === nothing + error("`SpeedMappingJL` requires `SpeedMapping.jl` to be loaded") + end + + return SpeedMappingJL(σ_min, stabilize, check_obj, orders, :SpeedMappingJL) +end + +""" + FixedPointAccelerationJL(; + algorithm = :Anderson, m = missing, condition_number_threshold = missing, + extrapolation_period = missing, replace_invalids = :NoAction + ) + +Wrapper over [FixedPointAcceleration.jl](https://s-baumann.github.io/FixedPointAcceleration.jl/) +for solving Fixed Point Problems. We allow using this algorithm to solve root finding +problems as well. + +### Keyword Arguments + + - `algorithm`: The algorithm to use. Can be `:Anderson`, `:MPE`, `:RRE`, `:VEA`, `:SEA`, + `:Simple`, `:Aitken` or `:Newton`. + - `m`: The number of previous iterates to use for the extrapolation. Only valid for + `:Anderson`. + - `condition_number_threshold`: The condition number threshold for Least Squares Problem. + Only valid for `:Anderson`. + - `extrapolation_period`: The number of iterates between extrapolations. Only valid for + `:MPE`, `:RRE`, `:VEA` and `:SEA`. Defaults to `7` for `:MPE` & `:RRE`, and `6` for + `:SEA` and `:VEA`. For `:SEA` and `:VEA`, this must be a multiple of `2`. + - `replace_invalids`: The method to use for replacing invalid iterates. Can be + `:ReplaceInvalids`, `:ReplaceVector` or `:NoAction`. + +!!! note + + This algorithm is only available if `FixedPointAcceleration.jl` is installed and loaded. +""" +@concrete struct FixedPointAccelerationJL <: AbstractNonlinearSolveAlgorithm + algorithm::Symbol + extrapolation_period::Int + replace_invalids::Symbol + dampening + m::Int + condition_number_threshold + name::Symbol +end + +function FixedPointAccelerationJL(; + algorithm = :Anderson, m = missing, condition_number_threshold = missing, + extrapolation_period = missing, replace_invalids = :NoAction, dampening = 1.0 +) + if Base.get_extension(@__MODULE__, :NonlinearSolveFixedPointAccelerationExt) === nothing + error("`FixedPointAccelerationJL` requires `FixedPointAcceleration.jl` to be loaded") + end + + @assert algorithm in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton) + @assert replace_invalids in (:ReplaceInvalids, :ReplaceVector, :NoAction) + + if algorithm !== :Anderson + @assert condition_number_threshold===missing "`condition_number_threshold` is only valid for Anderson acceleration" + @assert m===missing "`m` is only valid for Anderson acceleration" + end + condition_number_threshold === missing && (condition_number_threshold = 1e3) + m === missing && (m = 10) + + if algorithm !== :MPE && algorithm !== :RRE && algorithm !== :VEA && algorithm !== :SEA + @assert extrapolation_period===missing "`extrapolation_period` is only valid for MPE, RRE, VEA and SEA" + end + if extrapolation_period === missing + extrapolation_period = algorithm === :SEA || algorithm === :VEA ? 6 : 7 + else + if (algorithm === :SEA || algorithm === :VEA) && extrapolation_period % 2 != 0 + throw(AssertionError("`extrapolation_period` must be multiples of 2 for SEA and VEA")) + end + end + + return FixedPointAccelerationJL( + algorithm, extrapolation_period, replace_invalids, + dampening, m, condition_number_threshold, :FixedPointAccelerationJL + ) +end + +""" + SIAMFANLEquationsJL(; + method = :newton, delta = 1e-3, linsolve = nothing, autodiff = missing + ) + +### Keyword Arguments + + - `method`: the choice of method for solving the nonlinear system. + - `delta`: initial pseudo time step, default is 1e-3. + - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`. + - `m`: Depth for Anderson acceleration, default as 0 for Picard iteration. + - `beta`: Anderson mixing parameter, change f(x) to (1-beta)x+beta*f(x), + equivalent to accelerating damped Picard iteration. + - `autodiff`: Defaults to `missing`, which means we will default to letting + `SIAMFANLEquations` construct the jacobian if `f.jac` is not provided. In other cases, + we use it to generate a jacobian similar to other NonlinearSolve solvers. + +### Submethod Choice + + - `:newton`: Classical Newton method. + - `:pseudotransient`: Pseudo transient method. + - `:secant`: Secant method for scalar equations. + - `:anderson`: Anderson acceleration for fixed point iterations. + +!!! note + + This algorithm is only available if `SIAMFANLEquations.jl` is installed and loaded. +""" +@concrete struct SIAMFANLEquationsJL <: AbstractNonlinearSolveAlgorithm + method::Symbol + delta + linsolve <: Union{Symbol, Nothing} + m::Int + beta + autodiff + name::Symbol +end + +function SIAMFANLEquationsJL(; + method = :newton, delta = 1e-3, linsolve = nothing, m = 0, beta = 1.0, + autodiff = missing +) + if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing + error("`SIAMFANLEquationsJL` requires `SIAMFANLEquations.jl` to be loaded") + end + return SIAMFANLEquationsJL( + method, delta, linsolve, m, beta, autodiff, :SIAMFANLEquationsJL + ) +end + +""" + PETScSNES(; petsclib = missing, autodiff = nothing, mpi_comm = missing, kwargs...) + +Wrapper over [PETSc.jl](https://github.com/JuliaParallel/PETSc.jl) SNES solvers. + +### Keyword Arguments + + - `petsclib`: PETSc library to use. If set to `missing`, then we will use the first + available PETSc library in `PETSc.petsclibs` based on the problem element type. + - `autodiff`: the choice of method for generating the Jacobian. Defaults to `nothing` + which means that a default is selected according to the problem specification. Can be + any valid ADTypes.jl autodiff type (conditional on that backend being supported in + NonlinearSolve.jl). If set to `missing`, then PETSc computes the Jacobian using finite + differences. + - `mpi_comm`: MPI communicator to use. If set to `missing`, then we will use + `MPI.COMM_SELF`. + - `kwargs`: Keyword arguments to be passed to the PETSc SNES solver. See [PETSc SNES + documentation](https://petsc.org/release/manual/snes/) and + [SNESSetFromOptions](https://petsc.org/release/manualpages/SNES/SNESSetFromOptions) + for more information. + +### Options via `CommonSolve.solve` + +These options are forwarded from `solve` to the PETSc SNES solver. If these are provided to +`kwargs`, then they will be ignored. + +| `solve` option | PETSc SNES option | +|:-------------- |:----------------- | +| `atol` | `snes_atol` | +| `rtol` | `snes_rtol` | +| `maxiters` | `snes_max_it` | +| `show_trace` | `snes_monitor` | + +!!! note + + This algorithm is only available if `PETSc.jl` is installed and loaded. +""" +@concrete struct PETScSNES <: AbstractNonlinearSolveAlgorithm + petsclib + mpi_comm + autodiff + snes_options +end + +function PETScSNES(; petsclib = missing, autodiff = nothing, mpi_comm = missing, kwargs...) + if Base.get_extension(@__MODULE__, :NonlinearSolvePETScExt) === nothing + error("`PETScSNES` requires `PETSc.jl` to be loaded") + end + return PETScSNES(petsclib, mpi_comm, autodiff, kwargs) +end diff --git a/src/helpers.jl b/src/helpers.jl index e69de29bb..8b1378917 100644 --- a/src/helpers.jl +++ b/src/helpers.jl @@ -0,0 +1 @@ + diff --git a/src/internal/helpers.jl b/src/internal/helpers.jl deleted file mode 100644 index a7326bb51..000000000 --- a/src/internal/helpers.jl +++ /dev/null @@ -1,85 +0,0 @@ -# Extension Algorithm Helpers -function __test_termination_condition(termination_condition, alg) - !(termination_condition isa AbsNormTerminationMode) && - termination_condition !== nothing && - error("`$(alg)` does not support termination conditions!") -end - -function __construct_extension_f(prob::AbstractNonlinearProblem; alias_u0::Bool = false, - can_handle_oop::Val = False, can_handle_scalar::Val = False, - make_fixed_point::Val = False, force_oop::Val = False) - if can_handle_oop === False && can_handle_scalar === True - error("Incorrect Specification: OOP not supported but scalar supported.") - end - - resid = evaluate_f(prob, prob.u0) - u0 = can_handle_scalar === True || !(prob.u0 isa Number) ? - __maybe_unaliased(prob.u0, alias_u0) : [prob.u0] - - fₚ = if make_fixed_point === True - if isinplace(prob) - @closure (du, u) -> (prob.f(du, u, prob.p); du .+= u) - else - @closure u -> prob.f(u, prob.p) .+ u - end - else - if isinplace(prob) - @closure (du, u) -> prob.f(du, u, prob.p) - else - @closure u -> prob.f(u, prob.p) - end - end - - 𝐟 = if isinplace(prob) - u0_size, du_size = size(u0), size(resid) - @closure (du, u) -> (fₚ(reshape(du, du_size), reshape(u, u0_size)); du) - else - if prob.u0 isa Number - if can_handle_scalar === True - fₚ - elseif can_handle_oop === True - @closure u -> [fₚ(first(u))] - else - @closure (du, u) -> (du[1] = fₚ(first(u)); du) - end - else - u0_size = size(u0) - if can_handle_oop === True - @closure u -> vec(fₚ(reshape(u, u0_size))) - else - @closure (du, u) -> (copyto!(du, fₚ(reshape(u, u0_size))); du) - end - end - end - - 𝐅 = if force_oop === True && applicable(𝐟, u0, u0) - _resid = resid isa Number ? [resid] : _vec(resid) - du = _vec(zero(_resid)) - @closure u -> begin - 𝐟(du, u) - return du - end - else - 𝐟 - end - - return 𝐅, _vec(u0), (resid isa Number ? [resid] : _vec(resid)) -end - -function __construct_extension_jac(prob, alg, u0, fu; can_handle_oop::Val = False, - can_handle_scalar::Val = False, autodiff = nothing, initial_jacobian = False, - kwargs...) - autodiff = select_jacobian_autodiff(prob, autodiff) - - Jₚ = construct_jacobian_cache( - prob, alg, prob.f, fu, u0, prob.p; stats = empty_nlstats(), autodiff, kwargs...) - - 𝓙 = (can_handle_scalar === False && prob.u0 isa Number) ? @closure(u->[Jₚ(u[1])]) : Jₚ - - 𝐉 = (can_handle_oop === False && !isinplace(prob)) ? - @closure((J, u)->copyto!(J, 𝓙(u))) : 𝓙 - - initial_jacobian === False && return 𝐉 - - return 𝐉, Jₚ(nothing) -end diff --git a/test/23_test_problems_tests.jl b/test/23_test_problems_tests.jl index 6246cbecb..23b9a3e2a 100644 --- a/test/23_test_problems_tests.jl +++ b/test/23_test_problems_tests.jl @@ -1,144 +1,149 @@ -# @testsetup module RobustnessTesting -# using NonlinearSolve, LinearAlgebra, LinearSolve, NonlinearProblemLibrary, Test - -# problems = NonlinearProblemLibrary.problems -# dicts = NonlinearProblemLibrary.dicts - -# function test_on_library( -# problems, dicts, alg_ops, broken_tests, ϵ = 1e-4; skip_tests = nothing) -# for (idx, (problem, dict)) in enumerate(zip(problems, dicts)) -# x = dict["start"] -# res = similar(x) -# nlprob = NonlinearProblem(problem, copy(x)) -# @testset "$idx: $(dict["title"])" begin -# for alg in alg_ops -# try -# sol = solve(nlprob, alg; maxiters = 10000) -# problem(res, sol.u, nothing) - -# skip = skip_tests !== nothing && idx in skip_tests[alg] -# if skip -# @test_skip norm(res, Inf) ≤ ϵ -# continue -# end -# broken = idx in broken_tests[alg] ? true : false -# @test norm(res, Inf)≤ϵ broken=broken -# catch err -# @error err -# broken = idx in broken_tests[alg] ? true : false -# if broken -# @test false broken=true -# else -# @test 1 == 2 -# end -# end -# end -# end -# end -# end - -# export test_on_library, problems, dicts -# end - -# @testitem "PolyAlgorithms" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (RobustMultiNewton(), FastShortcutNonlinearPolyalg()) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [] -# broken_tests[alg_ops[2]] = [] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "NewtonRaphson" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (NewtonRaphson(),) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [1] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "TrustRegion" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple), -# TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan), -# TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei), -# TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan), -# TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin), -# TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.NLsolve)) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [11, 21] -# broken_tests[alg_ops[2]] = [11, 21] -# broken_tests[alg_ops[3]] = [11, 21] -# broken_tests[alg_ops[4]] = [8, 11, 21] -# broken_tests[alg_ops[5]] = [21] -# broken_tests[alg_ops[6]] = [11, 21] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "LevenbergMarquardt" setup=[RobustnessTesting] tags=[:core] begin -# using LinearSolve - -# alg_ops = (LevenbergMarquardt(), LevenbergMarquardt(; α_geodesic = 0.1), -# LevenbergMarquardt(; linsolve = CholeskyFactorization())) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [11, 21] -# broken_tests[alg_ops[2]] = [11, 21] -# broken_tests[alg_ops[3]] = [11, 21] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "DFSane" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (DFSane(),) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [1, 2, 3, 5, 21] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "Broyden" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)), -# Broyden(; update_rule = Val(:bad_broyden)), -# Broyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:bad_broyden))) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# if Sys.isapple() -# broken_tests[alg_ops[1]] = [1, 5, 11] -# broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18] -# broken_tests[alg_ops[3]] = [1, 5, 6, 9, 11] -# broken_tests[alg_ops[4]] = [5, 6, 8, 11] -# else -# broken_tests[alg_ops[1]] = [1, 5, 11, 15] -# broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18] -# broken_tests[alg_ops[3]] = [1, 5, 9, 11] -# broken_tests[alg_ops[4]] = [5, 6, 8, 11] -# end - -# test_on_library(problems, dicts, alg_ops, broken_tests, Sys.isapple() ? 1e-3 : 1e-4) -# end - -# @testitem "Klement" setup=[RobustnessTesting] tags=[:core] begin -# alg_ops = (Klement(), Klement(; init_jacobian = Val(:true_jacobian_diagonal))) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [1, 2, 4, 5, 11, 18, 22] -# broken_tests[alg_ops[2]] = [2, 4, 5, 7, 18, 22] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end - -# @testitem "PseudoTransient" setup=[RobustnessTesting] tags=[:core] begin -# # PT relies on the root being a stable equilibrium for convergence, so it won't work on -# # most problems -# alg_ops = (PseudoTransient(),) - -# broken_tests = Dict(alg => Int[] for alg in alg_ops) -# broken_tests[alg_ops[1]] = [1, 2, 3, 11, 15, 16] - -# test_on_library(problems, dicts, alg_ops, broken_tests) -# end +@testsetup module RobustnessTesting +using NonlinearSolve, LinearAlgebra, LinearSolve, NonlinearProblemLibrary, Test + +problems = NonlinearProblemLibrary.problems +dicts = NonlinearProblemLibrary.dicts + +function test_on_library( + problems, dicts, alg_ops, broken_tests, ϵ = 1e-4; skip_tests = nothing) + for (idx, (problem, dict)) in enumerate(zip(problems, dicts)) + x = dict["start"] + res = similar(x) + nlprob = NonlinearProblem(problem, copy(x)) + @testset "$idx: $(dict["title"])" begin + for alg in alg_ops + try + sol = solve(nlprob, alg; maxiters = 10000) + problem(res, sol.u, nothing) + + skip = skip_tests !== nothing && idx in skip_tests[alg] + if skip + @test_skip norm(res, Inf) ≤ ϵ + continue + end + broken = idx in broken_tests[alg] ? true : false + @test norm(res, Inf)≤ϵ broken=broken + catch err + @error err + broken = idx in broken_tests[alg] ? true : false + if broken + @test false broken=true + else + @test 1 == 2 + end + end + end + end + end +end + +export test_on_library, problems, dicts +end + +@testitem "PolyAlgorithms" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = (RobustMultiNewton(), FastShortcutNonlinearPolyalg()) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [] + broken_tests[alg_ops[2]] = [] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "NewtonRaphson" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = (NewtonRaphson(),) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [1] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "TrustRegion" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = ( + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Simple), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Fan), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Hei), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Yuan), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.Bastin), + TrustRegion(; radius_update_scheme = RadiusUpdateSchemes.NLsolve) + ) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [11, 21] + broken_tests[alg_ops[2]] = [11, 21] + broken_tests[alg_ops[3]] = [11, 21] + broken_tests[alg_ops[4]] = [8, 11, 21] + broken_tests[alg_ops[5]] = [21] + broken_tests[alg_ops[6]] = [11, 21] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "LevenbergMarquardt" setup=[RobustnessTesting] tags=[:core] begin + using LinearSolve + + alg_ops = ( + LevenbergMarquardt(), + LevenbergMarquardt(; α_geodesic = 0.1), + LevenbergMarquardt(; linsolve = CholeskyFactorization()) + ) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [11, 21] + broken_tests[alg_ops[2]] = [11, 21] + broken_tests[alg_ops[3]] = [11, 21] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "DFSane" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = (DFSane(),) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [1, 2, 3, 5, 21] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "Broyden" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = (Broyden(), Broyden(; init_jacobian = Val(:true_jacobian)), + Broyden(; update_rule = Val(:bad_broyden)), + Broyden(; init_jacobian = Val(:true_jacobian), update_rule = Val(:bad_broyden))) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + if Sys.isapple() + broken_tests[alg_ops[1]] = [1, 5, 11] + broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18] + broken_tests[alg_ops[3]] = [1, 5, 6, 9, 11] + broken_tests[alg_ops[4]] = [5, 6, 8, 11] + else + broken_tests[alg_ops[1]] = [1, 5, 11, 15] + broken_tests[alg_ops[2]] = [1, 5, 8, 11, 18] + broken_tests[alg_ops[3]] = [1, 5, 9, 11] + broken_tests[alg_ops[4]] = [5, 6, 8, 11] + end + + test_on_library(problems, dicts, alg_ops, broken_tests, Sys.isapple() ? 1e-3 : 1e-4) +end + +@testitem "Klement" setup=[RobustnessTesting] tags=[:core] begin + alg_ops = (Klement(), Klement(; init_jacobian = Val(:true_jacobian_diagonal))) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [1, 2, 4, 5, 11, 18, 22] + broken_tests[alg_ops[2]] = [2, 4, 5, 7, 18, 22] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end + +@testitem "PseudoTransient" setup=[RobustnessTesting] tags=[:core] begin + # PT relies on the root being a stable equilibrium for convergence, so it won't work on + # most problems + alg_ops = (PseudoTransient(),) + + broken_tests = Dict(alg => Int[] for alg in alg_ops) + broken_tests[alg_ops[1]] = [1, 2, 3, 11, 15, 16] + + test_on_library(problems, dicts, alg_ops, broken_tests) +end diff --git a/test/wrappers/fixedpoint_tests.jl b/test/wrappers/fixedpoint_tests.jl index a3ec497c7..a04e35069 100644 --- a/test/wrappers/fixedpoint_tests.jl +++ b/test/wrappers/fixedpoint_tests.jl @@ -1,11 +1,6 @@ -@testsetup module WrapperFixedPointImports -using Reexport -@reexport using LinearAlgebra -import SIAMFANLEquations, FixedPointAcceleration, SpeedMapping, NLsolve -end +@testitem "Simple Scalar Problem" tags=[:wrappers] begin + import SpeedMapping, SIAMFANLEquations, NLsolve, FixedPointAcceleration -# Simple Scalar Problem -@testitem "Simple Scalar Problem" setup=[WrapperFixedPointImports] tags=[:wrappers] begin f1(x, p) = cos(x) - x prob = NonlinearProblem(f1, 1.1) @@ -22,7 +17,9 @@ end end # Simple Vector Problem -@testitem "Simple Vector Problem" setup=[WrapperFixedPointImports] tags=[:wrappers] begin +@testitem "Simple Vector Problem" tags=[:wrappers] begin + import SpeedMapping, SIAMFANLEquations, NLsolve, FixedPointAcceleration + f2(x, p) = cos.(x) .- x prob = NonlinearProblem(f2, [1.1, 1.1]) @@ -41,7 +38,10 @@ end # Fixed Point for Power Method # Taken from https://github.com/NicolasL-S/SpeedMapping.jl/blob/95951db8f8a4457093090e18802ad382db1c76da/test/runtests.jl -@testitem "Power Method" setup=[WrapperFixedPointImports] tags=[:wrappers] begin +@testitem "Power Method" tags=[:wrappers] begin + using LinearAlgebra + import SpeedMapping, SIAMFANLEquations, NLsolve, FixedPointAcceleration + C = [1 2 3; 4 5 6; 7 8 9] A = C + C' B = Hermitian(ones(10) * ones(10)' .* im + Diagonal(1:10)) diff --git a/test/wrappers/least_squares_tests.jl b/test/wrappers/least_squares_tests.jl index 53cea758d..99d3e0ef3 100644 --- a/test/wrappers/least_squares_tests.jl +++ b/test/wrappers/least_squares_tests.jl @@ -1,52 +1,32 @@ @testsetup module WrapperNLLSSetup -using Reexport -@reexport using LinearAlgebra, StableRNGs, StaticArrays, Random, ForwardDiff, Zygote -import FastLevenbergMarquardt, LeastSquaresOptim, MINPACK -true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4]) -true_function(y, x, θ) = (@. y = θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])) +include("../../common/common_nlls_testing.jl") -θ_true = [1.0, 0.1, 2.0, 0.5] - -x = [-1.0, -0.5, 0.0, 0.5, 1.0] - -const y_target = true_function(x, θ_true) - -function loss_function(θ, p) - ŷ = true_function(p, θ) - return ŷ .- y_target -end - -function loss_function(resid, θ, p) - true_function(resid, p, θ) - resid .= resid .- y_target - return resid -end - -θ_init = θ_true .+ randn!(StableRNG(0), similar(θ_true)) * 0.1 - -export loss_function, θ_init, y_target, true_function, x, θ_true end @testitem "LeastSquaresOptim.jl" setup=[WrapperNLLSSetup] tags=[:wrappers] begin - prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x) - prob_iip = NonlinearLeastSquaresProblem( - NonlinearFunction(loss_function; resid_prototype = zero(y_target)), θ_init, x) + import LeastSquaresOptim nlls_problems = [prob_oop, prob_iip] - solvers = [LeastSquaresOptimJL(alg; autodiff) - for alg in (:lm, :dogleg), - autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff(), :central, :forward)] + solvers = [] + for alg in (:lm, :dogleg), + autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff(), :central, :forward) + + push!(solvers, LeastSquaresOptimJL(alg; autodiff)) + end for prob in nlls_problems, solver in solvers sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) @test SciMLBase.successful_retcode(sol) - @test norm(sol.resid, Inf) < 1e-6 + @test maximum(abs, sol.resid) < 1e-6 end end @testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Provided" setup=[WrapperNLLSSetup] tags=[:wrappers] begin + import FastLevenbergMarquardt, MINPACK + using ForwardDiff + function jac!(J, θ, p) resid = zeros(length(p)) ForwardDiff.jacobian!(J, (resid, θ) -> loss_function(resid, θ, p), resid, θ) @@ -58,19 +38,24 @@ end probs = [ NonlinearLeastSquaresProblem( NonlinearFunction{true}( - loss_function; resid_prototype = zero(y_target), jac = jac!), - θ_init, - x), + loss_function; resid_prototype = zero(y_target), jac = jac! + ), + θ_init, x + ), NonlinearLeastSquaresProblem( NonlinearFunction{false}( - loss_function; resid_prototype = zero(y_target), jac = jac), - θ_init, - x), + loss_function; resid_prototype = zero(y_target), jac = jac + ), + θ_init, x + ), NonlinearLeastSquaresProblem( - NonlinearFunction{false}(loss_function; jac), θ_init, x)] + NonlinearFunction{false}(loss_function; jac), θ_init, x + ) + ] solvers = Any[FastLevenbergMarquardtJL(linsolve) for linsolve in (:cholesky, :qr)] Sys.isapple() || push!(solvers, CMINPACK()) + for solver in solvers, prob in probs sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) @test maximum(abs, sol.resid) < 1e-6 @@ -78,28 +63,41 @@ end end @testitem "FastLevenbergMarquardt.jl + CMINPACK: Jacobian Not Provided" setup=[WrapperNLLSSetup] tags=[:wrappers] begin + import FastLevenbergMarquardt, MINPACK + probs = [ NonlinearLeastSquaresProblem( NonlinearFunction{true}(loss_function; resid_prototype = zero(y_target)), - θ_init, x), + θ_init, x + ), NonlinearLeastSquaresProblem( NonlinearFunction{false}(loss_function; resid_prototype = zero(y_target)), - θ_init, x), - NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function), θ_init, x)] + θ_init, x + ), + NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function), θ_init, x) + ] - solvers = vec(Any[FastLevenbergMarquardtJL(linsolve; autodiff) - for linsolve in (:cholesky, :qr), - autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff())]) - Sys.isapple() || - append!(solvers, [CMINPACK(; method) for method in (:auto, :lm, :lmdif)]) + solvers = [] + for linsolve in (:cholesky, :qr), + autodiff in (nothing, AutoForwardDiff(), AutoFiniteDiff()) + + push!(solvers, FastLevenbergMarquardtJL(linsolve; autodiff)) + end + if Sys.isapple() + for method in (:auto, :lm, :lmdif) + push!(solvers, CMINPACK(; method)) + end + end for solver in solvers, prob in probs sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8) - @test norm(sol.resid, Inf) < 1e-6 + @test maximum(abs, sol.resid) < 1e-6 end end @testitem "FastLevenbergMarquardt.jl + StaticArrays" setup=[WrapperNLLSSetup] tags=[:wrappers] begin + using StaticArrays, FastLevenbergMarquardt + x_sa = SA[-1.0, -0.5, 0.0, 0.5, 1.0] const y_target_sa = true_function(x_sa, θ_true) @@ -113,5 +111,5 @@ end prob_sa = NonlinearLeastSquaresProblem{false}(loss_function_sa, θ_init_sa, x) sol = solve(prob_sa, FastLevenbergMarquardtJL()) - @test norm(sol.resid, Inf) < 1e-6 + @test maximum(abs, sol.resid) < 1e-6 end diff --git a/test/wrappers/rootfind_tests.jl b/test/wrappers/rootfind_tests.jl index 852368ae3..4f317d76d 100644 --- a/test/wrappers/rootfind_tests.jl +++ b/test/wrappers/rootfind_tests.jl @@ -1,13 +1,6 @@ -@testsetup module WrapperRootfindImports -using Reexport -@reexport using LinearAlgebra -import NLSolvers, NLsolve, SIAMFANLEquations, MINPACK, PETSc +@testitem "Steady State Problems" tags=[:wrappers] begin + import NLSolvers, NLsolve, SIAMFANLEquations, MINPACK, PETSc -export NLSolvers -end - -@testitem "Steady State Problems" setup=[WrapperRootfindImports] tags=[:wrappers] begin - # IIP Tests function f_iip(du, u, p, t) du[1] = 2 - 2u[1] du[2] = u[1] - 4u[2] @@ -30,7 +23,6 @@ end @test maximum(abs, sol.resid) < 1e-6 end - # OOP Tests f_oop(u, p, t) = [2 - 2u[1], u[1] - 4u[2]] u0 = zeros(2) prob_oop = SteadyStateProblem(f_oop, u0) @@ -52,8 +44,10 @@ end end # Can lead to segfaults -@testitem "Nonlinear Root Finding Problems" setup=[WrapperRootfindImports] tags=[:wrappers] retries=3 begin - # IIP Tests +@testitem "Nonlinear Root Finding Problems" tags=[:wrappers] begin + using LinearAlgebra + import NLSolvers, NLsolve, SIAMFANLEquations, MINPACK, PETSc + function f_iip(du, u, p) du[1] = 2 - 2u[1] du[2] = u[1] - 4u[2] @@ -77,7 +71,6 @@ end @test maximum(abs, sol.resid) < 1e-6 end - # OOP Tests f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]] u0 = zeros(2) prob_oop = NonlinearProblem{false}(f_oop, u0) @@ -97,7 +90,6 @@ end @test maximum(abs, sol.resid) < 1e-6 end - # Tolerance Tests 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-15], @@ -156,9 +148,11 @@ end sol = solve(ProbN, NLsolveJL(); abstol = 1e-8) @test maximum(abs, sol.resid) < 1e-6 - sol = solve(ProbN, + sol = solve( + ProbN, NLSolversJL(NLSolvers.LineSearch(NLSolvers.Newton(), NLSolvers.Backtracking())); - abstol = 1e-8) + abstol = 1e-8 + ) @test maximum(abs, sol.resid) < 1e-6 sol = solve(ProbN, SIAMFANLEquationsJL(; method = :newton); abstol = 1e-8) @test maximum(abs, sol.resid) < 1e-6 @@ -170,7 +164,9 @@ end end end -@testitem "PETSc SNES Floating Points" setup=[WrapperRootfindImports] tags=[:wrappers] skip=:(Sys.iswindows()) begin +@testitem "PETSc SNES Floating Points" tags=[:wrappers] skip=:(Sys.iswindows()) begin + import PETSc + f(u, p) = u .* u .- 2 u0 = [1.0, 1.0]