10x - 100x speed up when using own implementation of Implicit Euler or Implicit Midpoint on linear test problem #2586
Comments
|
I assume this is because |
|
If you check |
|
Note that there is also a significant difference with iterative solvers (code copied from the OP): # using not the default linear solvers my_implicit_euler still 10x faster
solver = KrylovJL_GMRES();
@btime sol4 = my_implicit_euler(prob_op, dt, solver);
@btime sol5 = solve(prob_op, ImplicitEuler(linsolve = solver), dt=dt, adaptive = false, save_everystep = false);
@btime sol6 = solve(prob, ImplicitEuler(linsolve = solver), dt=dt, adaptive = false, save_everystep = false); |
julia> sol2.stats
SciMLBase.DEStats
Number of function 1 evaluations: 510005
Number of function 2 evaluations: 0
Number of W matrix evaluations: 5000
Number of linear solves: 5004
Number of Jacobians created: 5000
Number of nonlinear solver iterations: 5004
Number of nonlinear solver convergence failures: 0
Number of fixed-point solver iterations: 0
Number of fixed-point solver convergence failures: 0
Number of rootfind condition calls: 0
Number of accepted steps: 5000
Number of rejected steps: 0 |
|
Oh adaptive=false doesn't have a mechanism for rejections and storing factorizations. |
|
This is now slowly beyond me expertise, but does this also explain the 10x difference when using indirect solvers like |
|
Maybe I am wrong, but isn't the W matrix only needed when solving nonlinear systems, and could it be that considerable time is spent evaluating it even when solving linear problems? julia> sol6 = solve(prob, ImplicitEuler(linsolve = KrylovJL_GMRES()), dt=0.01, adaptive = false, save_everystep = false);
julia> sol6.stats
SciMLBase.DEStats
Number of function 1 evaluations: 1507
Number of function 2 evaluations: 0
Number of W matrix evaluations: 500
Number of linear solves: 506
Number of Jacobians created: 0
Number of nonlinear solver iterations: 506
Number of nonlinear solver convergence failures: 0
Number of fixed-point solver iterations: 0
Number of fixed-point solver convergence failures: 0
Number of rootfind condition calls: 0
Number of accepted steps: 500
Number of rejected steps: 0Or maybe it is just different default tolerances in |
|
You also need it for linear systems, but for a linear ODE you have that W is always the same, and so lu(W) is always the same. What's going on here is that W is constructed 500 times, and lu(W) is done 500 times, when you only need to do it once. Again, this is a biproduct of safety measures on adaptive=false, and we can improve that. |
|
But there is also a big difference when no LU decomposition is used, just an iterative solver. |
|
Can you share a profile for that case? |
|
This is the result of running # profiling
solver = KrylovJL_GMRES();
# my implicit euler
@profview for i in 1:500 my_implicit_euler(prob_op, dt, solver); end
# ImplicitEuler
@profview for i in 1:500 solve(prob_op, ImplicitEuler(linsolve = solver), dt=dt, adaptive = false, save_everystep = false); endhttps://github.com/cwittens/two_html_profiles I couldn't upload .html files to this comment and pictures or just text makes to profiles not so nice to look at I figured, so I uploaded them to a repo. |
|
it's a difference in which krylov solver we're using. We're spending all the time here https://github.com/JuliaSmoothOptimizers/Krylov.jl/blob/f03bafbf3d7124fb7397493741a7420f39b82042/src/krylov_solvers.jl#L2500C5-L2500C15 which seems very odd. We appear to be making thousands of tiny arrays instead of a matrix or something like that... |
|
Note sure if this is relevant/related, but the same performance different can be seen using And ┌ Warning: cg! doesn't support right preconditioning.
└ @ LinearSolve C:\Users\colli\.julia\packages\LinearSolve\WDeMC\src\iterative_wrappers.jl:286
┌ Warning: cg! doesn't support right preconditioning.
└ @ LinearSolve C:\Users\colli\.julia\packages\LinearSolve\WDeMC\src\iterative_wrappers.jl:286
┌ Warning: cg! doesn't support right preconditioning.
└ @ LinearSolve C:\Users\colli\.julia\packages\LinearSolve\WDeMC\src\iterative_wrappers.jl:286
....even when explicitly telling the FunctionOperator that it is posdef and symmetric |
Describe the bug 🐞
When using my own implementation of Implicit Euler or Implicit Midpoint (definition below), I get a 10x - 100x speed up on the linear test problem
(du, u, p, t) -> du .= p .* uwhen usingFunctionOperatorin the ODEProblem definition. (and also when solving Parabolic PDEs as discovered by @ranocha , but this is not shown in this MRE).Minimal Reproducible Example 👇
Output:
And the results are pretty much the same for my_implicit_midpoint and ImplicitMidpoint.
Output:
Obviously the speed up depends on the system size.
Environment (please complete the following information):
using Pkg; Pkg.status()using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)versioninfo()Additional context
Add any other context about the problem here.
The text was updated successfully, but these errors were encountered: