-
-
Notifications
You must be signed in to change notification settings - Fork 210
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Change to LinearSolve.jl #1543
Change to LinearSolve.jl #1543
Conversation
Test case is: using OrdinaryDiffEq, StaticArrays, BenchmarkTools
function rober(du,u,p,t)
y₁,y₂,y₃ = u
k₁,k₂,k₃ = p
@inbounds begin
du[1] = -k₁*y₁+k₃*y₂*y₃
du[2] = k₁*y₁-k₂*y₂^2-k₃*y₂*y₃
du[3] = k₂*y₂^2
end
nothing
end
prob = ODEProblem(rober,[1.0,0.0,0.0],(0.0,1e5),[0.04,3e7,1e4])
@btime solve(prob,Rosenbrock23()) # 761.578 μs (633 allocations: 53.05 KiB) Only the Rosenbrock and Extrapolation methods will be tedious, the rest should already be handled by the Newton part. |
src/OrdinaryDiffEq.jl
Outdated
@@ -176,6 +178,7 @@ using DocStringExtensions | |||
include("interp_func.jl") | |||
include("composite_algs.jl") | |||
|
|||
#= |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
Need to not forget to undo this 😅
src/caches/rosenbrock_caches.jl
Outdated
@@ -73,7 +73,9 @@ function alg_cache(alg::Rosenbrock23,u,rate_prototype,::Type{uEltypeNoUnits},::T | |||
tf = TimeGradientWrapper(f,uprev,p) | |||
uf = UJacobianWrapper(f,t,p) | |||
linsolve_tmp = zero(rate_prototype) | |||
linsolve = alg.linsolve(Val{:init},uf,u) | |||
|
|||
linprob = LinearProblem(W,vec(u); u0=vec(u)) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
need to do this to the whole file, and the whole extrapolation caches file
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
sounds cumbersome. vim macros are might speed it up
src/nlsolve/newton.jl
Outdated
end | ||
linsolve = set_b(linsolve,b) | ||
linres = solve(linsolve,weights=weight,reltol=reltol) | ||
copyto!(dz,linres.u) |
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
This adds one copyto!
operation, but it shouldn't be bad.
It's the first time I see this error message. Does it still occur? |
fc96127
to
0266d59
Compare
Currently hits a weird error: ```julia [ Info: Precompiling OrdinaryDiffEq [1dea7af3-3e70-54e6-95c3-0bf5283fa5ed] ERROR: LoadError: syntax: SSAValue objects should not occur in an AST Stacktrace: [1] top-level scope @ C:\Users\accou\.julia\dev\OrdinaryDiffEq\src\perform_step\rosenbrock_perform_step.jl:27 [2] include(mod::Module, _path::String) @ Base .\Base.jl:418 [3] include(x::String) @ OrdinaryDiffEq C:\Users\accou\.julia\dev\OrdinaryDiffEq\src\OrdinaryDiffEq.jl:4 [4] top-level scope @ C:\Users\accou\.julia\dev\OrdinaryDiffEq\src\OrdinaryDiffEq.jl:148 [5] include @ .\Base.jl:418 [inlined] [6] include_package_for_output(pkg::Base.PkgId, input::String, depot_path::Vector{String}, dl_load_path::Vector{String}, load_path::Vector{String}, concrete_deps::Vector{Pair{Base.PkgId, UInt64}}, source::String) @ Base .\loading.jl:1318 [7] top-level scope @ none:1 [8] eval @ .\boot.jl:373 [inlined] [9] eval(x::Expr) @ Base.MainInclude .\client.jl:453 [10] top-level scope @ none:1 in expression starting at C:\Users\accou\.julia\dev\OrdinaryDiffEq\src\perform_step\rosenbrock_perform_step.jl:27 in expression starting at C:\Users\accou\.julia\dev\OrdinaryDiffEq\src\OrdinaryDiffEq.jl:1 ``` @YingboMa do you know what could cause that?
Tests are passing, so now I'm just going to refactor a bit to condense the code before merging 🎉 |
Benchmarks look to be comfortably a wash even though it allocates more. using OrdinaryDiffEq, SnoopCompile, ForwardDiff
lorenz = (du,u,p,t) -> begin
du[1] = 10.0(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]; tspan = (0.0,100.0);
prob = ODEProblem(lorenz,u0,tspan); alg = Rodas5();
tinf = @snoopi_deep ForwardDiff.gradient(u0 -> sum(solve(ODEProblem(lorenz,u0,tspan),alg)), u0)
tinf = @snoopi_deep ForwardDiff.gradient(u0 -> sum(solve(ODEProblem(lorenz,u0,tspan),alg)), u0)
# Before
InferenceTimingNode: 1.219960/3.583559 on Core.Compiler.Timings.ROOT() with 34 direct children
InferenceTimingNode: 1.064226/1.756057 on Core.Compiler.Timings.ROOT() with 13 direct children
# After
InferenceTimingNode: 1.257616/3.774263 on Core.Compiler.Timings.ROOT() with 33 direct children
InferenceTimingNode: 1.006660/1.756162 on Core.Compiler.Timings.ROOT() with 13 direct children
using OrdinaryDiffEq, SnoopCompile, ForwardDiff, BenchmarkTools
lorenz = (du,u,p,t) -> begin
du[1] = 10.0(u[2]-u[1])
du[2] = u[1]*(28.0-u[3]) - u[2]
du[3] = u[1]*u[2] - (8/3)*u[3]
end
u0 = [1.0;0.0;0.0]; tspan = (0.0,100.0);
prob = ODEProblem(lorenz,u0,tspan); alg = Rodas5();
tinf = @snoopi_deep solve(ODEProblem(lorenz,u0,tspan),alg)
tinf = @snoopi_deep solve(ODEProblem(lorenz,u0,tspan),alg)
@btime solve(ODEProblem(lorenz,u0,tspan),alg)
@btime solve(ODEProblem(lorenz,u0,(0.0,10.0)),alg)
# Before
InferenceTimingNode: 1.060946/2.763417 on Core.Compiler.Timings.ROOT() with 8 direct children
InferenceTimingNode: 0.002333/0.002333 on Core.Compiler.Timings.ROOT() with 0 direct children
2.055 ms (14069 allocations: 781.14 KiB)
175.000 μs (1341 allocations: 75.97 KiB)
# After
InferenceTimingNode: 1.129433/3.072010 on Core.Compiler.Timings.ROOT() with 7 direct children
InferenceTimingNode: 0.002182/0.002182 on Core.Compiler.Timings.ROOT() with 0 direct children
2.221 ms (22155 allocations: 2.24 MiB)
162.900 μs (1917 allocations: 184.00 KiB) ROBER Before after Pollu Before after |
SciML/DelayDiffEq.jl#222 Downstream PRs are already ready. |
Currently hits a weird error:
@YingboMa do you know what could cause that?