Merge pull request #1 from JuliaComputing/myb/benchmark

Modify the benchmark and tests

Author: YingboMa

src/NonlinearSolve.jl

@@ -27,7 +27,6 @@ module NonlinearSolve
 
   # DiffEq styled algorithms
   export Bisection, Falsi, NewtonRaphson
-  export ScalarBisection, ScalarNewton
 
   export reinit!
 end # module

src/scalar.jl

@@ -1,11 +1,4 @@
-"""
-ScalarNewton
-
-Fast Newton Raphson for scalar problems.
-"""
-struct ScalarNewton <: AbstractNonlinearSolveAlgorithm end
-
-function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarNewton, args...; xatol = nothing, xrtol = nothing, maxiters = 1000, kwargs...) where {uType}
+function DiffEqBase.solve(prob::NonlinearProblem{<:Number}, ::NewtonRaphson, args...; xatol = nothing, xrtol = nothing, maxiters = 1000, kwargs...)
   f = Base.Fix2(prob.f, prob.p)
   x = float(prob.u0)
   T = typeof(x)
@@ -15,26 +8,18 @@ function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarNewton,
   xo = oftype(x, Inf)
   for i in 1:maxiters
     fx, dfx = value_derivative(f, x)
-    iszero(fx) && return x
+    iszero(fx) && return NewtonSolution(x, :Default)
     Δx = dfx \ fx
     x -= Δx
     if isapprox(x, xo, atol=atol, rtol=rtol)
-        return x
+        return NewtonSolution(x, :Default)
     end
     xo = x
   end
-  return oftype(x, NaN)
+  return NewtonSolution(x, :MaxitersExceeded)
 end
 
-"""
-  ScalarBisection
-
-Fast Bisection for scalar problems. Note that it doesn't returns exact solution, but returns
-the best left limit of the exact solution.
-"""
-struct ScalarBisection <: AbstractNonlinearSolveAlgorithm end
-
-function DiffEqBase.solve(prob::NonlinearProblem{uType, false}, ::ScalarBisection, args...; maxiters = 1000, kwargs...) where {uType}
+function DiffEqBase.solve(prob::NonlinearProblem, ::Bisection, args...; maxiters = 1000, kwargs...)
   f = Base.Fix2(prob.f, prob.p)
   left, right = prob.u0
   fl, fr = f(left), f(right)

src/solve.jl

@@ -9,10 +9,11 @@ end
 function DiffEqBase.init(prob::NonlinearProblem{uType, iip}, alg::AbstractBracketingAlgorithm, args...;
     alias_u0 = false,
     maxiters = 1000,
+    # bracketing algorithms only solve scalar problems
     immutable = (prob.u0 isa StaticArray || prob.u0 isa Number),
     kwargs...
   ) where {uType, iip}
-  
+
   if !(prob.u0 isa Tuple)
     error("You need to pass a tuple of u0 in bracketing algorithms.")
   end
@@ -60,7 +61,7 @@ function DiffEqBase.init(prob::NonlinearProblem{uType, iip}, alg::AbstractNewton
     fu = f(u, p)
   end
 
-  
+
   sol = build_newton_solution(u, Val(iip))
   if immutable
     return NewtonImmutableSolver(1, f, alg, u, fu, p, nothing, false, maxiters, internalnorm, :Default, tol, sol)
@@ -81,7 +82,7 @@ function DiffEqBase.solve!(solver::AbstractNonlinearSolver)
   if solver.iter == solver.maxiters
     solver.retcode = :MaxitersExceeded
   end
-  set_solution!(solver)
+  solver = set_solution(solver)
   return solver.sol
 end
 
@@ -151,19 +152,25 @@ function check_for_exact_solution!(solver::BracketingSolver)
   return false
 end
 
-function set_solution!(solver::BracketingSolver)
-  solver.sol.left = solver.left
-  solver.sol.right = solver.right
-  solver.sol.retcode = solver.retcode
+function set_solution(solver::BracketingSolver)
+  sol = solver.sol
+  @set! sol.left = solver.left
+  @set! sol.right = solver.right
+  @set! sol.retcode = solver.retcode
+  @set! solver.sol = sol
+  return solver
 end
 
 function get_solution(solver::BracketingImmutableSolver)
   return (left = solver.left, right = solver.right, retcode = solver.retcode)
 end
 
-function set_solution!(solver::NewtonSolver)
-  solver.sol.u = solver.u
-  solver.sol.retcode = solver.retcode
+function set_solution(solver::NewtonSolver)
+  sol = solver.sol
+  @set! sol.u = solver.u
+  @set! sol.retcode = solver.retcode
+  @set! solver.sol = sol
+  return solver
 end
 
 function get_solution(solver::NewtonImmutableSolver)

src/types.jl

@@ -83,7 +83,7 @@ function build_solution(u_prototype, ::Val{false})
     return BracketingSolution(zero(u_prototype), zero(u_prototype), :Default)
 end
 
-mutable struct NewtonSolution{uType}
+struct NewtonSolution{uType}
     u::uType
     retcode::Symbol
 end

test/runtests.jl

@@ -3,26 +3,35 @@ using StaticArrays
 using BenchmarkTools
 using Test
 
-function benchmark_immutable()
-    probN = NonlinearProblem((u,p) -> u .* u .- 2, @SVector[1.0, 1.0])
+function benchmark_immutable(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
     solver = init(probN, NewtonRaphson(), immutable = true, tol = 1e-9)
-    sol = @btime solve!($solver)
-    @test all(sol.u .* sol.u .- 2 .< 1e-9)
+    sol = solve!(solver)
 end
 
-function benchmark_mutable()
-    probN = NonlinearProblem((u,p) -> u .* u .- 2, @SVector[1.0, 1.0])
+function benchmark_mutable(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
     solver = init(probN, NewtonRaphson(), immutable = false, tol = 1e-9)
-    sol = @btime (reinit!($solver, $probN); solve!($solver))
-    @test all(sol.u .* sol.u .- 2 .< 1e-9)
+    sol = (reinit!(solver, probN); solve!(solver))
 end
 
-function benchmark_scalar()
-    probN = NonlinearProblem((u,p) -> u .* u .- 2, 1.0)
-    sol = @btime (solve($probN, ScalarNewton()))
-    @test sol * sol - 2 < 1e-9
+function benchmark_scalar(f, u0)
+    probN = NonlinearProblem{false}(f, u0)
+    sol = (solve(probN, NewtonRaphson()))
 end
 
-benchmark_immutable()
-benchmark_mutable()
-benchmark_scalar()
+f, u0 = (u,p) -> u .* u .- 2, @SVector[1.0, 1.0]
+sf, su0 = (u,p) -> u * u - 2, 1.0
+sol = benchmark_immutable(f, u0)
+@test sol.retcode === :Default
+@test all(sol.u .* sol.u .- 2 .< 1e-9)
+sol = benchmark_mutable(f, u0)
+@test sol.retcode === :Default
+@test all(sol.u .* sol.u .- 2 .< 1e-9)
+sol = benchmark_scalar(sf, su0)
+@test sol.retcode === :Default
+@test sol.u * sol.u - 2 < 1e-9
+
+@test (@ballocated benchmark_immutable($f, $u0)) == 0
+@test (@ballocated benchmark_mutable($f, $u0)) < 200
+@test (@ballocated benchmark_scalar($sf, $su0)) == 0