diff --git a/.github/workflows/Downstream.yml b/.github/workflows/Downstream.yml index 916a33b4b..274675fcd 100644 --- a/.github/workflows/Downstream.yml +++ b/.github/workflows/Downstream.yml @@ -18,6 +18,8 @@ jobs: os: [ubuntu-latest] package: - {user: SciML, repo: ModelingToolkit.jl, group: All} + - {user: SciML, repo: OrdinaryDiffEq.jl, group: Interface} + - {user: SciML, repo: OrdinaryDiffEq.jl, group: Regression} steps: - uses: actions/checkout@v2 - uses: julia-actions/setup-julia@v1 diff --git a/Project.toml b/Project.toml index b7bd102f6..79c2f12ec 100644 --- a/Project.toml +++ b/Project.toml @@ -1,17 +1,20 @@ name = "NonlinearSolve" uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec" -authors = ["Kanav Gupta "] +authors = ["SciML"] version = "0.3.24" [deps] ArrayInterfaceCore = "30b0a656-2188-435a-8636-2ec0e6a096e2" +DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e" FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae" RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd" Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" +SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" @@ -19,9 +22,10 @@ UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed" ArrayInterfaceCore = "0.1.1" FiniteDiff = "2" ForwardDiff = "0.10.3" +LinearSolve = "1" RecursiveArrayTools = "2" Reexport = "0.2, 1" -SciMLBase = "1.32" +SciMLBase = "1.73" Setfield = "0.7, 0.8, 1" StaticArrays = "0.12,1.0" UnPack = "1.0" diff --git a/docs/src/basics/FAQ.md b/docs/src/basics/FAQ.md index 9c521e57d..0038ed50b 100644 --- a/docs/src/basics/FAQ.md +++ b/docs/src/basics/FAQ.md @@ -15,15 +15,15 @@ myfun(x, lv) = x * sin(x) - lv function f(out, levels, u0) for i in 1:N - out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun), + out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), u0, levels[i]), Falsi()).u end end function f2(out, levels, u0) for i in 1:N - out[i] = solve(NonlinearProblem{false}(NonlinearFunction{false}(myfun), - u0, levels[i]), NewtonRaphson()).u + out[i] = solve(IntervalNonlinearProblem{false}(IntervalNonlinearFunction{false}(myfun), + u0, levels[i]), SimpleNewtonRaphson()).u end end diff --git a/docs/src/basics/NonlinearFunctions.md b/docs/src/basics/NonlinearFunctions.md index b420d891b..c02520e2c 100644 --- a/docs/src/basics/NonlinearFunctions.md +++ b/docs/src/basics/NonlinearFunctions.md @@ -9,5 +9,6 @@ of pre-computed functions to speed up the calculations. This is offered via the ## Function Type Definitions ```@docs +SciMLBase.IntervalNonlinearFunction SciMLBase.NonlinearFunction ``` diff --git a/docs/src/basics/NonlinearProblem.md b/docs/src/basics/NonlinearProblem.md index dda50c5c7..33cfd7b80 100644 --- a/docs/src/basics/NonlinearProblem.md +++ b/docs/src/basics/NonlinearProblem.md @@ -1,5 +1,32 @@ # Nonlinear Problems +## The Two Types of Nonlinear Problems + +NonlinearSolve.jl tackles two related types of nonlinear systems: + +1. Interval rootfinding problems. I.e., find the ``t in [t_0, t_f]`` such that `f(t) = 0`. +2. Systems of nonlinear equations, i.e. find the `u` such that `f(u) = 0`. + +The former is for solving scalar rootfinding problems, i.e. finding a single number, and +requires that a bracketing interval is known. For a bracketing interval, one must have that +the sign of `f(t_0)` is opposite the sign of `f(t_f)`, thus guaranteeing a root in the +interval. + +The latter type of nonlinear system can be multidimensional and thus no ordering nor +boundaries are assumed to be known. For a system of nonlinear equations, `f` can return +an array and the solver seeks to find the value of `u` for which all outputs of `f` are +simultaniously zero. + +!!! note + + Interval rootfinding problems allow for `f` to return an array, in which case the interval + rootfinding problem is interpreted as finding the first `t` such that any of the components + of the array hit zero. + + +## Problem Construction Details + ```@docs +SciMLBase.IntervalNonlinearProblem SciMLBase.NonlinearProblem ``` diff --git a/docs/src/index.md b/docs/src/index.md index eea1a76f8..35419a3c2 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -12,8 +12,7 @@ allow for automatically generating high-performance code. Performance is key: the current methods are made to be highly performant on scalar and statically sized small problems, with options for large-scale systems. -If you run into any performance issues, please file an issue. Note that this -package is distinct from [SciMLNLSolve.jl](https://github.com/SciML/SciMLNLSolve.jl). +If you run into any performance issues, please file an issue. Consult the [NonlinearSystemSolvers](@ref nonlinearsystemsolvers) page for information on how to import solvers from different packages. @@ -79,7 +78,7 @@ Pkg.status(;mode = PKGMODE_MANIFEST) # hide ``` ```@raw html -You can also download the +You can also download the project file. -``` \ No newline at end of file +``` diff --git a/docs/src/solvers/BracketingSolvers.md b/docs/src/solvers/BracketingSolvers.md index f83650029..0dd6c4375 100644 --- a/docs/src/solvers/BracketingSolvers.md +++ b/docs/src/solvers/BracketingSolvers.md @@ -1,8 +1,8 @@ -# Bracketing Solvers +# Interval Rootfinding Methods (Bracketing Solvers) -`solve(prob::NonlinearProblem,alg;kwargs)` +`solve(prob::IntervalNonlinearProblem,alg;kwargs)` -Solves for ``f(u)=0`` in the problem defined by `prob` using the algorithm +Solves for ``f(t)=0`` in the problem defined by `prob` using the algorithm `alg`. If no algorithm is given, a default algorithm will be chosen. This page is solely focused on the bracketing methods for scalar nonlinear equations. @@ -14,7 +14,10 @@ less stable than `Bisection`. ## Full List of Methods -### NonlinearSolve.jl +### SimpleNonlinearSolve.jl + +These methods are automatically included as part of NonlinearSolve.jl. Though one can use +SimpleNonlinearSolve.jl directly to decrease the dependencies and improve load time. - `Falsi`: A non-allocating regula falsi method - `Bisection`: A common bisection method diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md index 6e9206a59..7602c1936 100644 --- a/docs/src/solvers/NonlinearSystemSolvers.md +++ b/docs/src/solvers/NonlinearSystemSolvers.md @@ -9,13 +9,15 @@ This page is solely focused on the methods for nonlinear systems. ## Recommended Methods -`NewtonRaphson` is a good choice for most problems. It is non-allocating on -static arrays and thus really well-optimized for small systems, while for large +`NewtonRaphson` is a good choice for most problems. For large systems it can make use of sparsity patterns for sparse automatic differentiation and sparse linear solving of very large systems. That said, as a classic Newton -method, its stability region can be smaller than other methods. `NLSolveJL`'s +method, its stability region can be smaller than other methods. Meanwhile, `SimpleNewtonRaphson` +is an implementation which is specialized for small equations. It is non-allocating on +static arrays and thus really well-optimized for small systems, thus usually outperforming +the other methods when such types are used for `u0`. `NLSolveJL`'s `:trust_region` method can be a good choice for high stability, along with -`CMINPACK`. +`CMINPACK`.s For a system which is very non-stiff (i.e., the condition number of the Jacobian is small, or the eigenvalues of the Jacobian are within a few orders of magnitude), @@ -29,8 +31,17 @@ These are the core solvers. - `NewtonRaphson(;autodiff=true,chunk_size=12,diff_type=Val{:forward},linsolve=DEFAULT_LINSOLVE)`: A Newton-Raphson method with swappable nonlinear solvers and autodiff methods - for high performance on large and sparse systems. When used on objects like - static arrays, this method is non-allocating. + for high performance on large and sparse systems. + +### SimpleNonlinearSolve.jl + +These methods are included with NonlinearSolve.jl by default, though SimpleNonlinearSolve.jl +can be used directly to reduce dependencies and improve load times. + +- `SimpleNewtonRaphson()`: A simplified implementation of the Newton-Raphson method. Has the + property that when used with states `u` as a `Number` or `StaticArray`, the solver is + very efficient and non-allocating. Thus this implmentation is well-suited for small + systems of equations. ### SciMLNLSolve.jl diff --git a/docs/src/tutorials/iterator_interface.md b/docs/src/tutorials/iterator_interface.md index f66c6c768..9e1e70898 100644 --- a/docs/src/tutorials/iterator_interface.md +++ b/docs/src/tutorials/iterator_interface.md @@ -4,11 +4,8 @@ There is an iterator form of the nonlinear solver which mirrors the DiffEq integ ```julia f(u, p) = u .* u .- 2.0 -u0 = (1.0, 2.0) # brackets +u0 = 1.5 # brackets probB = NonlinearProblem(f, u0) -solver = init(probB, Falsi()) # Can iterate the solver object -solver = solve!(solver) +cache = init(probB, NewtonRaphson()) # Can iterate the solver object +solver = solve!(cache) ``` - -Note that the `solver` object is actually immutable since we want to make it -live on the stack for the sake of performance. diff --git a/docs/src/tutorials/nonlinear.md b/docs/src/tutorials/nonlinear.md index 73007f441..45be3ffb1 100644 --- a/docs/src/tutorials/nonlinear.md +++ b/docs/src/tutorials/nonlinear.md @@ -27,6 +27,6 @@ a bracket instead of an initial condition, for example: ```julia f(u, p) = u .* u .- 2.0 u0 = (1.0, 2.0) # brackets -probB = NonlinearProblem(f, u0) +probB = IntervalNonlinearProblem(f, u0) sol = solve(probB, Falsi()) ``` diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl index a1fbd9b4e..df53d4bd5 100644 --- a/src/NonlinearSolve.jl +++ b/src/NonlinearSolve.jl @@ -9,25 +9,28 @@ using StaticArrays using RecursiveArrayTools using LinearAlgebra import ArrayInterfaceCore +import LinearSolve +using DiffEqBase @reexport using SciMLBase +@reexport using SimpleNonlinearSolve abstract type AbstractNonlinearSolveAlgorithm <: SciMLBase.AbstractNonlinearAlgorithm end -abstract type AbstractBracketingAlgorithm <: AbstractNonlinearSolveAlgorithm end abstract type AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ} <: AbstractNonlinearSolveAlgorithm end -abstract type AbstractImmutableNonlinearSolver <: AbstractNonlinearSolveAlgorithm end + +function SciMLBase.__solve(prob::NonlinearProblem, + alg::AbstractNonlinearSolveAlgorithm, args...; + kwargs...) + cache = init(prob, alg, args...; kwargs...) + sol = solve!(cache) +end include("utils.jl") include("jacobian.jl") -include("types.jl") -include("solve.jl") -include("bisection.jl") -include("falsi.jl") include("raphson.jl") -include("scalar.jl") +include("ad.jl") -# DiffEq styled algorithms -export Bisection, Falsi, NewtonRaphson +export NewtonRaphson end # module diff --git a/src/ad.jl b/src/ad.jl new file mode 100644 index 000000000..1e310f9f9 --- /dev/null +++ b/src/ad.jl @@ -0,0 +1,39 @@ +function scalar_nlsolve_ad(prob, alg, args...; kwargs...) + f = prob.f + p = value(prob.p) + + u0 = value(prob.u0) + newprob = NonlinearProblem(f, u0, p; prob.kwargs...) + + sol = solve(newprob, alg, args...; kwargs...) + + uu = sol.u + if p isa Number + f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p) + else + f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p) + end + + f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu) + pp = prob.p + sumfun = let f_x′ = -f_x + ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p) + end + partials = sum(sumfun, zip(f_p, pp)) + return sol, partials +end + +function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip, + <:Dual{T, V, P}}, alg::NewtonRaphson, + args...; kwargs...) where {iip, T, V, P} + sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) + return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; + retcode = sol.retcode) +end +function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip, + <:AbstractArray{<:Dual{T, V, P}}}, + alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P} + sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) + return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; + retcode = sol.retcode) +end diff --git a/src/bisection.jl b/src/bisection.jl deleted file mode 100644 index 7268da9f1..000000000 --- a/src/bisection.jl +++ /dev/null @@ -1,116 +0,0 @@ -struct Bisection <: AbstractBracketingAlgorithm - exact_left::Bool - exact_right::Bool -end - -function Bisection(; exact_left = false, exact_right = false) - Bisection(exact_left, exact_right) -end - -struct BisectionCache{uType} - state::Int - left::uType - right::uType -end - -function alg_cache(alg::Bisection, left, right, p, ::Val{true}) - BisectionCache(0, left, right) -end - -function alg_cache(alg::Bisection, left, right, p, ::Val{false}) - BisectionCache(0, left, right) -end - -function perform_step(solver::BracketingImmutableSolver, alg::Bisection, cache) - @unpack f, p, left, right, fl, fr, cache = solver - - if cache.state == 0 - fzero = zero(fl) - fl * fr > fzero && - error("Bracket became non-containing in between iterations. This could mean that " - * - "input function crosses the x axis multiple times. Bisection is not the right method to solve this.") - - mid = (left + right) / 2 - - if left == mid || right == mid - @set! solver.force_stop = true - @set! solver.retcode = ReturnCode.FloatingPointLimit - return solver - end - - fm = f(mid, p) - - if iszero(fm) - if alg.exact_left - @set! cache.state = 1 - @set! cache.right = mid - @set! cache.left = mid - @set! solver.cache = cache - elseif alg.exact_right - @set! solver.left = prevfloat_tdir(mid, left, right) - solver = sync_residuals!(solver) - @set! cache.state = 2 - @set! cache.left = mid - @set! solver.cache = cache - else - @set! solver.left = prevfloat_tdir(mid, left, right) - @set! solver.right = nextfloat_tdir(mid, left, right) - solver = sync_residuals!(solver) - @set! solver.force_stop = true - return solver - end - else - if sign(fm) == sign(fl) - @set! solver.left = mid - @set! solver.fl = fm - else - @set! solver.right = mid - @set! solver.fr = fm - end - end - elseif cache.state == 1 - mid = (left + cache.right) / 2 - - if cache.right == mid || left == mid - if alg.exact_right - @set! cache.state = 2 - @set! solver.cache = cache - return solver - else - @set! solver.right = nextfloat_tdir(mid, left, right) - solver = sync_residuals!(solver) - @set! solver.force_stop = true - return solver - end - end - - fm = f(mid, p) - - if iszero(fm) - @set! cache.right = mid - @set! solver.cache = cache - else - @set! solver.left = mid - @set! solver.fl = fm - end - else - mid = (cache.left + right) / 2 - - if right == mid || cache.left == mid - @set! solver.force_stop = true - return solver - end - - fm = f(mid, p) - - if iszero(fm) - @set! cache.left = mid - @set! solver.cache = cache - else - @set! solver.right = mid - @set! solver.fr = fm - end - end - solver -end diff --git a/src/falsi.jl b/src/falsi.jl deleted file mode 100644 index f3420cc63..000000000 --- a/src/falsi.jl +++ /dev/null @@ -1,46 +0,0 @@ -struct Falsi <: AbstractBracketingAlgorithm end - -function alg_cache(alg::Falsi, left, right, p, ::Val{true}) - nothing -end - -function alg_cache(alg::Falsi, left, right, p, ::Val{false}) - nothing -end - -function perform_step(solver, alg::Falsi, cache) - @unpack f, p, left, right, fl, fr = solver - - fzero = zero(fl) - fl * fr > fzero && - error("Bracket became non-containing in between iterations. This could mean that " - * - "input function crosses the x axis multiple times. Bisection is not the right method to solve this.") - - mid = (fr * left - fl * right) / (fr - fl) - - if right == mid || right == mid - @set! solver.force_stop = true - @set! solver.retcode = ReturnCode.FloatingPointLimit - return solver - end - - fm = f(mid, p) - - if iszero(fm) - # todo: phase 2 bisection similar to the raw method - @set! solver.force_stop = true - @set! solver.left = mid - @set! solver.fl = fm - @set! solver.retcode = ReturnCode.ExactSolutionLeft - else - if sign(fm) == sign(fl) - @set! solver.left = mid - @set! solver.fl = fm - else - @set! solver.right = mid - @set! solver.fr = fm - end - end - return solver -end diff --git a/src/jacobian.jl b/src/jacobian.jl index f10f61137..96c07ada9 100644 --- a/src/jacobian.jl +++ b/src/jacobian.jl @@ -28,11 +28,21 @@ function calc_J(solver, uf::ImmutableJacobianWrapper) return J end +function jacobian(f, x::Number, solver) + if alg_autodiff(solver.alg) + J = ForwardDiff.derivative(f, x) + else + J = FiniteDiff.finite_difference_derivative(f, x, alg_difftype(solver.alg), + eltype(x)) + end + return J +end + function jacobian(f, x, solver) if alg_autodiff(solver.alg) J = ForwardDiff.jacobian(f, x) else - J = FiniteDiff.finite_difference_jacobian(f, x, solver.alg.diff_type, eltype(x)) + J = FiniteDiff.finite_difference_jacobian(f, x, alg_difftype(solver.alg), eltype(x)) end return J end diff --git a/src/raphson.jl b/src/raphson.jl index 4866e61ab..b1362596d 100644 --- a/src/raphson.jl +++ b/src/raphson.jl @@ -4,8 +4,6 @@ struct NewtonRaphson{CS, AD, FDT, L, P, ST, CJ} <: precs::P end -DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing - function NewtonRaphson(; chunk_size = Val{0}(), autodiff = Val{true}(), standardtag = Val{true}(), concrete_jac = nothing, diff_type = Val{:forward}, linsolve = nothing, precs = DEFAULT_PRECS) @@ -14,68 +12,55 @@ function NewtonRaphson(; chunk_size = Val{0}(), autodiff = Val{true}(), _unwrap_val(concrete_jac)}(linsolve, precs) end -mutable struct NewtonRaphsonCache{ufType, L, jType, uType, JC} +mutable struct NewtonRaphsonCache{iip, fType, algType, uType, duType, resType, pType, + INType, tolType, + probType, ufType, L, jType, JC} + f::fType + alg::algType + u::uType + fu::resType + p::pType uf::ufType linsolve::L J::jType - du1::uType + du1::duType jac_config::JC -end - -function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing, - du = nothing, u = nothing, p = nothing, t = nothing, - weight = nothing, solverdata = nothing, - reltol = nothing) where {P} - A !== nothing && (linsolve = LinearSolve.set_A(linsolve, A)) - b !== nothing && (linsolve = LinearSolve.set_b(linsolve, b)) - linu !== nothing && (linsolve = LinearSolve.set_u(linsolve, linu)) - - Plprev = linsolve.Pl isa LinearSolve.ComposePreconditioner ? linsolve.Pl.outer : - linsolve.Pl - Prprev = linsolve.Pr isa LinearSolve.ComposePreconditioner ? linsolve.Pr.outer : - linsolve.Pr - - _Pl, _Pr = precs(linsolve.A, du, u, p, nothing, A !== nothing, Plprev, Prprev, - solverdata) - if (_Pl !== nothing || _Pr !== nothing) - _weight = weight === nothing ? - (linsolve.Pr isa Diagonal ? linsolve.Pr.diag : linsolve.Pr.inner.diag) : - weight - Pl, Pr = wrapprecs(_Pl, _Pr, _weight) - linsolve = LinearSolve.set_prec(linsolve, Pl, Pr) - end - - linres = if reltol === nothing - solve(linsolve; reltol) - else - solve(linsolve; reltol) - end - - return linres -end - -function wrapprecs(_Pl, _Pr, weight) - if _Pl !== nothing - Pl = LinearSolve.ComposePreconditioner(LinearSolve.InvPreconditioner(Diagonal(_vec(weight))), - _Pl) - else - Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))) - end - - if _Pr !== nothing - Pr = LinearSolve.ComposePreconditioner(Diagonal(_vec(weight)), _Pr) - else - Pr = Diagonal(_vec(weight)) + iter::Int + force_stop::Bool + maxiters::Int + internalnorm::INType + retcode::SciMLBase.ReturnCode.T + abstol::tolType + prob::probType + + function NewtonRaphsonCache{iip}(f::fType, alg::algType, u::uType, fu::resType, + p::pType, + uf::ufType, linsolve::L, J::jType, du1::duType, + jac_config::JC, iter::Int, + force_stop::Bool, maxiters::Int, internalnorm::INType, + retcode::SciMLBase.ReturnCode.T, abstol::tolType, + prob::probType) where { + iip, fType, algType, uType, + duType, resType, pType, INType, + tolType, + probType, ufType, L, jType, JC} + new{iip, fType, algType, uType, duType, resType, pType, INType, tolType, + probType, ufType, L, jType, JC}(f, alg, u, fu, p, + uf, linsolve, J, du1, jac_config, iter, + force_stop, maxiters, internalnorm, + retcode, abstol, prob) end - Pl, Pr end -function alg_cache(alg::NewtonRaphson, f, u, p, ::Val{true}) +function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{true}) uf = JacobianWrapper(f, p) J = false .* u .* u' - linprob = LinearProblem(W, _vec(zero(u)); u0 = _vec(zero(u))) - Pl, Pr = wrapprecs(alg.precs(W, nothing, u, p, nothing, nothing, nothing, nothing, + linprob = LinearProblem(J, _vec(zero(u)); u0 = _vec(zero(u))) + weight = similar(u) + recursivefill!(weight, false) + + Pl, Pr = wrapprecs(alg.precs(J, nothing, u, p, nothing, nothing, nothing, nothing, nothing)..., weight) linsolve = init(linprob, alg.linsolve, alias_A = true, alias_b = true, Pl = Pl, Pr = Pr) @@ -94,38 +79,79 @@ function alg_cache(alg::NewtonRaphson, f, u, p, ::Val{true}) alg.diff_type, eltype(u)) end end - NewtonRaphsonCache(uf, linsolve, J, du1, jac_config) + uf, linsolve, J, du1, jac_config +end + +function jacobian_caches(alg::NewtonRaphson, f, u, p, ::Val{false}) + nothing, nothing, nothing, nothing, nothing end -function alg_cache(alg::NewtonRaphson, f, u, p, ::Val{false}) - nothing +function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::NewtonRaphson, + args...; + alias_u0 = false, + maxiters = 1000, + abstol = 1e-6, + internalnorm = DEFAULT_NORM, + kwargs...) where {uType, iip} + if alias_u0 + u = prob.u0 + else + u = deepcopy(prob.u0) + end + f = prob.f + p = prob.p + if iip + fu = zero(u) + f(fu, u, p) + else + fu = f(u, p) + end + uf, linsolve, J, du1, jac_config = jacobian_caches(alg, f, u, p, Val(iip)) + + return NewtonRaphsonCache{iip}(f, alg, u, fu, p, uf, linsolve, J, du1, jac_config, + 1, false, maxiters, internalnorm, + ReturnCode.Default, abstol, prob) end -function perform_step(solver::NewtonImmutableSolver, alg::NewtonRaphson, ::Val{true}) - @unpack u, fu, f, p, cache = solver +function perform_step!(cache::NewtonRaphsonCache{true}) + @unpack u, fu, f, p, cache = cache @unpack J, linsolve, du1 = cache - calc_J!(J, solver, cache) + calc_J!(J, cache, cache) # u = u - J \ fu - linsolve = dolinsolve(alg.precs, solver.linsolve, A = J, b = fu, u = du1, - p = p, reltol = solver.tol) - @set! cache.linsolve = linsolve + linsolve = dolinsolve(alg.precs, linsolve, A = J, b = fu, u = du1, + p = p, reltol = cache.abstol) + cache.linsolve = linsolve @. u = u - du1 f(fu, u, p) - if solver.internalnorm(solver.fu) < solver.tol - @set! solver.force_stop = true + if cache.internalnorm(cache.fu) < cache.abstol + cache.force_stop = true + end + return nothing +end + +function perform_step!(cache::NewtonRaphsonCache{false}) + @unpack u, fu, f, p = cache + J = calc_J(cache, ImmutableJacobianWrapper(f, p)) + cache.u = u - J \ fu + fu = f(cache.u, p) + cache.fu = fu + if iszero(cache.fu) || cache.internalnorm(cache.fu) < cache.abstol + cache.force_stop = true end - return solver + return nothing end -function perform_step(solver::NewtonImmutableSolver, alg::NewtonRaphson, ::Val{false}) - @unpack u, fu, f, p = solver - J = calc_J(solver, ImmutableJacobianWrapper(f, p)) - @set! solver.u = u - J \ fu - fu = f(solver.u, p) - @set! solver.fu = fu - if iszero(solver.fu) || solver.internalnorm(solver.fu) < solver.tol - @set! solver.force_stop = true +function SciMLBase.solve!(cache::NewtonRaphsonCache) + while !cache.force_stop && cache.iter < cache.maxiters + perform_step!(cache) + cache.iter += 1 end - return solver + + if cache.iter == cache.maxiters + cache.retcode = ReturnCode.MaxIters + end + + SciMLBase.build_solution(cache.prob, cache.alg, cache.u, cache.fu; + retcode = cache.retcode) end diff --git a/src/scalar.jl b/src/scalar.jl deleted file mode 100644 index 07f2c392b..000000000 --- a/src/scalar.jl +++ /dev/null @@ -1,226 +0,0 @@ -function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}}, - alg::NewtonRaphson, args...; xatol = nothing, xrtol = nothing, - maxiters = 1000, kwargs...) - f = Base.Fix2(prob.f, prob.p) - x = float(prob.u0) - fx = float(prob.u0) - T = typeof(x) - atol = xatol !== nothing ? xatol : oneunit(eltype(T)) * (eps(one(eltype(T))))^(4 // 5) - rtol = xrtol !== nothing ? xrtol : eps(one(eltype(T)))^(4 // 5) - - if typeof(x) <: Number - xo = oftype(one(eltype(x)), Inf) - else - xo = map(x -> oftype(one(eltype(x)), Inf), x) - end - - for i in 1:maxiters - if alg_autodiff(alg) - fx, dfx = value_derivative(f, x) - elseif x isa AbstractArray - fx = f(x) - dfx = FiniteDiff.finite_difference_jacobian(f, x, alg.diff_type, eltype(x), fx) - else - fx = f(x) - dfx = FiniteDiff.finite_difference_derivative(f, x, alg.diff_type, eltype(x), - fx) - end - iszero(fx) && - return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Default) - Δx = dfx \ fx - x -= Δx - if isapprox(x, xo, atol = atol, rtol = rtol) - return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Default) - end - xo = x - end - return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters) -end - -function scalar_nlsolve_ad(prob, alg, args...; kwargs...) - f = prob.f - p = value(prob.p) - u0 = value(prob.u0) - - newprob = NonlinearProblem(f, u0, p; prob.kwargs...) - sol = solve(newprob, alg, args...; kwargs...) - - uu = sol.u - if p isa Number - f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p) - else - f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p) - end - - f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu) - pp = prob.p - sumfun = let f_x′ = -f_x - ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p) - end - partials = sum(sumfun, zip(f_p, pp)) - return sol, partials -end - -function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip, - <:Dual{T, V, P}}, alg::NewtonRaphson, - args...; kwargs...) where {iip, T, V, P} - sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) - return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) -end -function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector}, iip, - <:AbstractArray{<:Dual{T, V, P}}}, - alg::NewtonRaphson, args...; kwargs...) where {iip, T, V, P} - sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) - return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid; - retcode = sol.retcode) -end - -# avoid ambiguities -for Alg in [Bisection] - @eval function SciMLBase.solve(prob::NonlinearProblem{uType, iip, <:Dual{T, V, P}}, - alg::$Alg, args...; - kwargs...) where {uType, iip, T, V, P} - sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) - return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), - sol.resid; retcode = sol.retcode, - left = Dual{T, V, P}(sol.left, partials), - right = Dual{T, V, P}(sol.right, partials)) - #return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode, sol.resid) - end - @eval function SciMLBase.solve(prob::NonlinearProblem{uType, iip, - <:AbstractArray{<:Dual{T, V, P}}}, - alg::$Alg, args...; - kwargs...) where {uType, iip, T, V, P} - sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...) - return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), - sol.resid; retcode = sol.retcode, - left = Dual{T, V, P}(sol.left, partials), - right = Dual{T, V, P}(sol.right, partials)) - #return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode, sol.resid) - end -end - -function SciMLBase.solve(prob::NonlinearProblem, alg::Bisection, args...; maxiters = 1000, - kwargs...) - f = Base.Fix2(prob.f, prob.p) - left, right = prob.u0 - fl, fr = f(left), f(right) - - if iszero(fl) - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.ExactSolutionLeft, left = left, - right = right) - end - - i = 1 - if !iszero(fr) - while i < maxiters - mid = (left + right) / 2 - (mid == left || mid == right) && - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.FloatingPointLimit, - left = left, right = right) - fm = f(mid) - if iszero(fm) - right = mid - break - end - if sign(fl) == sign(fm) - fl = fm - left = mid - else - fr = fm - right = mid - end - i += 1 - end - end - - while i < maxiters - mid = (left + right) / 2 - (mid == left || mid == right) && - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.FloatingPointLimit, - left = left, right = right) - fm = f(mid) - if iszero(fm) - right = mid - fr = fm - else - left = mid - fl = fm - end - i += 1 - end - - return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters, - left = left, right = right) -end - -function SciMLBase.solve(prob::NonlinearProblem, alg::Falsi, args...; maxiters = 1000, - kwargs...) - f = Base.Fix2(prob.f, prob.p) - left, right = prob.u0 - fl, fr = f(left), f(right) - - if iszero(fl) - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.ExactSolutionLeft, left = left, - right = right) - end - - i = 1 - if !iszero(fr) - while i < maxiters - if nextfloat_tdir(left, prob.u0...) == right - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.FloatingPointLimit, - left = left, right = right) - end - mid = (fr * left - fl * right) / (fr - fl) - for i in 1:10 - mid = max_tdir(left, prevfloat_tdir(mid, prob.u0...), prob.u0...) - end - if mid == right || mid == left - break - end - fm = f(mid) - if iszero(fm) - right = mid - break - end - if sign(fl) == sign(fm) - fl = fm - left = mid - else - fr = fm - right = mid - end - i += 1 - end - end - - while i < maxiters - mid = (left + right) / 2 - (mid == left || mid == right) && - return SciMLBase.build_solution(prob, alg, left, fl; - retcode = ReturnCode.FloatingPointLimit, - left = left, right = right) - fm = f(mid) - if iszero(fm) - right = mid - fr = fm - elseif sign(fm) == sign(fl) - left = mid - fl = fm - else - right = mid - fr = fm - end - i += 1 - end - - return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters, - left = left, right = right) -end diff --git a/src/solve.jl b/src/solve.jl deleted file mode 100644 index 930ce9c7e..000000000 --- a/src/solve.jl +++ /dev/null @@ -1,123 +0,0 @@ -function SciMLBase.__solve(prob::NonlinearProblem, - alg::AbstractNonlinearSolveAlgorithm, args...; - kwargs...) - solver = init(prob, alg, args...; kwargs...) - sol = solve!(solver) -end - -function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, - alg::AbstractBracketingAlgorithm, args...; - alias_u0 = false, - maxiters = 1000, - kwargs...) where {uType, iip} - if !(prob.u0 isa Tuple) - error("You need to pass a tuple of u0 in bracketing algorithms.") - end - - if eltype(prob.u0) isa AbstractArray - error("Bracketing Algorithms work for scalar arguments only") - end - - if alias_u0 - left, right = prob.u0 - else - left, right = deepcopy(prob.u0) - end - f = prob.f - p = prob.p - fl = f(left, p) - fr = f(right, p) - cache = alg_cache(alg, left, right, p, Val(iip)) - return BracketingImmutableSolver(1, f, alg, left, right, fl, fr, p, false, maxiters, - ReturnCode.Default, cache, iip, prob) -end - -function SciMLBase.__init(prob::NonlinearProblem{uType, iip}, alg::AbstractNewtonAlgorithm, - args...; - alias_u0 = false, - maxiters = 1000, - tol = 1e-6, - internalnorm = DEFAULT_NORM, - kwargs...) where {uType, iip} - if alias_u0 - u = prob.u0 - else - u = deepcopy(prob.u0) - end - f = prob.f - p = prob.p - if iip - fu = zero(u) - f(fu, u, p) - else - fu = f(u, p) - end - cache = alg_cache(alg, f, u, p, Val(iip)) - return NewtonImmutableSolver(1, f, alg, u, fu, p, false, maxiters, internalnorm, - Retcode.Default, tol, cache, iip, prob) -end - -function SciMLBase.solve!(solver::AbstractImmutableNonlinearSolver) - solver = mic_check(solver) - while !solver.force_stop && solver.iter < solver.maxiters - solver = perform_step(solver, solver.alg, Val(solver.iip)) - @set! solver.iter += 1 - end - if solver.iter == solver.maxiters - @set! solver.retcode = ReturnCode.MaxIters - end - if typeof(solver) <: NewtonImmutableSolver - SciMLBase.build_solution(solver.prob, solver.alg, solver.u, solver.fu; - retcode = solver.retcode) - else - SciMLBase.build_solution(solver.prob, solver.alg, solver.left, solver.fl; - retcode = solver.retcode, left = solver.left, - right = solver.right) - end -end - -""" - mic_check(solver::AbstractImmutableNonlinearSolver) - mic_check!(solver::AbstractNonlinearSolver) - -Checks before running main solving iterations. -""" -function mic_check(solver::BracketingImmutableSolver) - @unpack f, fl, fr = solver - flr = fl * fr - fzero = zero(flr) - (flr > fzero) && error("Non bracketing interval passed in bracketing method.") - if fl == fzero - @set! solver.force_stop = true - @set! solver.retcode = Retcode.ExactSolutionLeft - elseif fr == fzero - @set! solver.force_stop = true - @set! solver.retcode = Retcode.ExactionSolutionRight - end - solver -end - -function mic_check(solver::NewtonImmutableSolver) - solver -end - -""" - reinit!(solver, prob) - -Reinitialize solver to the original starting conditions -""" -function SciMLBase.reinit!(solver::NewtonImmutableSolver, - prob::NonlinearProblem{uType, true}) where {uType} - @. solver.u = prob.u0 - @set! solver.iter = 1 - @set! solver.force_stop = false - return solver -end - -function SciMLBase.reinit!(solver::NewtonImmutableSolver, - prob::NonlinearProblem{uType, false}) where {uType} - @set! solver.u = prob.u0 - @set! solver.iter = 1 - @set! solver.force_stop = false - return solver -end diff --git a/src/types.jl b/src/types.jl deleted file mode 100644 index 99ecadba2..000000000 --- a/src/types.jl +++ /dev/null @@ -1,51 +0,0 @@ -struct BracketingImmutableSolver{fType, algType, uType, resType, pType, cacheType, probType - } <: AbstractImmutableNonlinearSolver - iter::Int - f::fType - alg::algType - left::uType - right::uType - fl::resType - fr::resType - p::pType - force_stop::Bool - maxiters::Int - retcode::SciMLBase.ReturnCode.T - cache::cacheType - iip::Bool - prob::probType -end - -# function BracketingImmutableSolver(iip, iter, f, alg, left, right, fl, fr, p, force_stop, maxiters, retcode, cache) -# BracketingImmutableSolver{iip, typeof(f), typeof(alg), -# typeof(left), typeof(fl), typeof(p), typeof(cache)}(iter, f, alg, left, right, fl, fr, p, force_stop, maxiters, retcode, cache) -# end - -struct NewtonImmutableSolver{fType, algType, uType, resType, pType, INType, tolType, - cacheType, probType} <: AbstractImmutableNonlinearSolver - iter::Int - f::fType - alg::algType - u::uType - fu::resType - p::pType - force_stop::Bool - maxiters::Int - internalnorm::INType - retcode::SciMLBase.ReturnCode.T - tol::tolType - cache::cacheType - iip::Bool - prob::probType -end - -# function NewtonImmutableSolver{iip}(iter, f, alg, u, fu, p, force_stop, maxiters, internalnorm, retcode, tol, cache) where iip -# NewtonImmutableSolver{iip, typeof(f), typeof(alg), typeof(u), -# typeof(fu), typeof(p), typeof(internalnorm), typeof(tol), typeof(cache)}(iter, f, alg, u, fu, p, force_stop, maxiters, internalnorm, retcode, tol, cache) -# end - -function sync_residuals!(solver::BracketingImmutableSolver) - @set! solver.fl = solver.f(solver.left, solver.p) - @set! solver.fr = solver.f(solver.right, solver.p) - solver -end diff --git a/src/utils.jl b/src/utils.jl index cf61b9afc..318ca01f9 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -151,3 +151,65 @@ value_derivative(f::F, x::SVector) where {F} = f(x), ForwardDiff.jacobian(f, x) value(x) = x value(x::Dual) = ForwardDiff.value(x) value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x) + +_unwrap_val(::Val{B}) where {B} = B +_unwrap_val(B) = B + +_vec(v) = vec(v) +_vec(v::Number) = v +_vec(v::AbstractVector) = v + +function alg_difftype(alg::AbstractNewtonAlgorithm{CS, AD, FDT, ST, CJ}) where {CS, AD, FDT, + ST, CJ} + FDT +end + +DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, cachedata) = nothing, nothing + +function dolinsolve(precs::P, linsolve; A = nothing, linu = nothing, b = nothing, + du = nothing, u = nothing, p = nothing, t = nothing, + weight = nothing, cachedata = nothing, + reltol = nothing) where {P} + A !== nothing && (linsolve = LinearSolve.set_A(linsolve, A)) + b !== nothing && (linsolve = LinearSolve.set_b(linsolve, b)) + linu !== nothing && (linsolve = LinearSolve.set_u(linsolve, linu)) + + Plprev = linsolve.Pl isa LinearSolve.ComposePreconditioner ? linsolve.Pl.outer : + linsolve.Pl + Prprev = linsolve.Pr isa LinearSolve.ComposePreconditioner ? linsolve.Pr.outer : + linsolve.Pr + + _Pl, _Pr = precs(linsolve.A, du, u, p, nothing, A !== nothing, Plprev, Prprev, + cachedata) + if (_Pl !== nothing || _Pr !== nothing) + _weight = weight === nothing ? + (linsolve.Pr isa Diagonal ? linsolve.Pr.diag : linsolve.Pr.inner.diag) : + weight + Pl, Pr = wrapprecs(_Pl, _Pr, _weight) + linsolve = LinearSolve.set_prec(linsolve, Pl, Pr) + end + + linres = if reltol === nothing + solve(linsolve) + else + solve(linsolve; reltol) + end + + return linres +end + +function wrapprecs(_Pl, _Pr, weight) + if _Pl !== nothing + Pl = LinearSolve.ComposePreconditioner(LinearSolve.InvPreconditioner(Diagonal(_vec(weight))), + _Pl) + else + Pl = LinearSolve.InvPreconditioner(Diagonal(_vec(weight))) + end + + if _Pr !== nothing + Pr = LinearSolve.ComposePreconditioner(Diagonal(_vec(weight)), _Pr) + else + Pr = Diagonal(_vec(weight)) + end + Pl, Pr +end diff --git a/test/basictests.jl b/test/basictests.jl index 98ee8627b..e4ab88103 100644 --- a/test/basictests.jl +++ b/test/basictests.jl @@ -5,14 +5,14 @@ using Test function benchmark_immutable(f, u0) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, NewtonRaphson(), tol = 1e-9) + solver = init(probN, NewtonRaphson(), abstol = 1e-9) sol = solve!(solver) end function benchmark_mutable(f, u0) probN = NonlinearProblem{false}(f, u0) - solver = init(probN, NewtonRaphson(), tol = 1e-9) - sol = (reinit!(solver, probN); solve!(solver)) + solver = init(probN, NewtonRaphson(), abstol = 1e-9) + sol = solve!(solver) end function benchmark_scalar(f, u0) @@ -39,9 +39,9 @@ sol = benchmark_scalar(sf, csu0) @test sol.retcode === ReturnCode.Default @test sol.u * sol.u - 2 < 1e-9 -@test (@ballocated benchmark_immutable(ff, cu0)) == 0 +@test (@ballocated benchmark_immutable(ff, cu0)) < 200 @test (@ballocated benchmark_mutable(ff, cu0)) < 200 -@test (@ballocated benchmark_scalar(sf, csu0)) == 0 +@test (@ballocated benchmark_scalar(sf, csu0)) < 400 # AD Tests using ForwardDiff @@ -51,7 +51,7 @@ f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0] g = function (p) probN = NonlinearProblem{false}(f, csu0, p) - sol = solve(probN, NewtonRaphson(), tol = 1e-9) + sol = solve(probN, NewtonRaphson(), abstol = 1e-9) return sol.u[end] end @@ -66,7 +66,7 @@ f, u0 = (u, p) -> u * u - p, 1.0 # NewtonRaphson g = function (p) probN = NonlinearProblem{false}(f, oftype(p, u0), p) - sol = solve(probN, NewtonRaphson()) + sol = solve(probN, NewtonRaphson(), abstol = 1e-10) return sol.u end @@ -77,34 +77,9 @@ for p in 1.1:0.1:100.0 @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) end -u0 = (1.0, 20.0) -# Falsi -g = function (p) - probN = NonlinearProblem{false}(f, typeof(p).(u0), p) - sol = solve(probN, Falsi()) - return sol.left -end - -for p in 1.1:0.1:100.0 - @test g(p) ≈ sqrt(p) - @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p)) -end - -f, u0 = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0) +f = (u, p) -> p[1] * u * u - p[2] t = (p) -> [sqrt(p[2] / p[1])] p = [0.9, 50.0] -for alg in [Bisection(), Falsi()] - global g, p - g = function (p) - probN = NonlinearProblem{false}(f, u0, p) - sol = solve(probN, Bisection()) - return [sol.left] - end - - @test g(p) ≈ [sqrt(p[2] / p[1])] - @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p) -end - gnewton = function (p) probN = NonlinearProblem{false}(f, 0.5, p) sol = solve(probN, NewtonRaphson()) @@ -119,8 +94,6 @@ f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0] probN = NonlinearProblem(f, u0) @test solve(probN, NewtonRaphson()).u[end] ≈ sqrt(2.0) -@test solve(probN, NewtonRaphson(); immutable = false).u[end] ≈ sqrt(2.0) -@test solve(probN, NewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) @test solve(probN, NewtonRaphson(; autodiff = false)).u[end] ≈ sqrt(2.0) for u0 in [1.0, [1, 1.0]] @@ -133,49 +106,3 @@ for u0 in [1.0, [1, 1.0]] @test solve(probN, NewtonRaphson()).u ≈ sol @test solve(probN, NewtonRaphson(; autodiff = false)).u ≈ sol end - -# Bisection Tests -f, u0 = (u, p) -> u .* u .- 2.0, (1.0, 2.0) -probB = NonlinearProblem(f, u0) - -# Falsi -solver = init(probB, Falsi()) -sol = solve!(solver) -@test sol.left ≈ sqrt(2.0) - -# this should call the fast scalar overload -@test solve(probB, Bisection()).left ≈ sqrt(2.0) - -# these should call the iterator version -solver = init(probB, Bisection()) -@test solver isa NonlinearSolve.BracketingImmutableSolver -@test solve!(solver).left ≈ sqrt(2.0) - -# Garuntee Tests for Bisection -f = function (u, p) - if u < 2.0 - return u - 2.0 - elseif u > 3.0 - return u - 3.0 - else - return 0.0 - end -end -probB = NonlinearProblem(f, (0.0, 4.0)) - -solver = init(probB, Bisection(; exact_left = true)) -sol = solve!(solver) -@test f(sol.left, nothing) < 0.0 -@test f(nextfloat(sol.left), nothing) >= 0.0 - -solver = init(probB, Bisection(; exact_right = true)) -sol = solve!(solver) -@test f(sol.right, nothing) > 0.0 -@test f(prevfloat(sol.right), nothing) <= 0.0 - -solver = init(probB, Bisection(; exact_left = true, exact_right = true); immutable = false) -sol = solve!(solver) -@test f(sol.left, nothing) < 0.0 -@test f(nextfloat(sol.left), nothing) >= 0.0 -@test f(sol.right, nothing) > 0.0 -@test f(prevfloat(sol.right), nothing) <= 0.0 diff --git a/test/runtests.jl b/test/runtests.jl index 0cf5ed83b..7269ed6bf 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -7,7 +7,6 @@ const is_APPVEYOR = Sys.iswindows() && haskey(ENV, "APPVEYOR") @time begin -if GROUP == "All" || GROUP == "Interface" - #@time @safetestset "Linear Solver Tests" begin include("interface/linear_solver_test.jl") end +if GROUP == "All" || GROUP == "Core" @time @safetestset "Basic Tests + Some AD" begin include("basictests.jl") end end end