Would you be willing to update the code in the benchmark to
condition1(u, t, integrator) = (t-0.5) * (t-0.5005)
cb1 = ContinuousCallback(condition1, nothing; rootfind=true,save_positions=(true,false), abstol=0.001)
const Mask2 = Ref(false)
function condition2(u, t, integrator)
if Mask2[] == false
if u[1] < 0
-1.
else
u[2]
end
else
-1.
end
end
function affect2(integrator)
Mask2[] = true
end
cb2 = ContinuousCallback(condition2, affect2; rootfind=true,save_positions=(true,false),abstol=1e-9)
const Mask3 = Ref(false)
function condition3(u, t, integrator)
if Mask3[] == false
if u[1] > 0
-1.
else
u[2]
end
else
-1.
end
end
function affect3(integrator)
integrator.u = X0
Mask3[] = true
end
and update the README to reflect the conclusions? (2x difference to same speed depending on the julia solver you use). When benchmarking vs competitors, I think it is important to make fair comparisons. (Imo, using the StaticArray version would be better, but thatâs a less big deal).
Yes, I will update the script to use static arrays as this definitely improves the performance of Julia code. Rest of your changes hardly have any effect or even negative as I showed in the numbers earlier. Of course, we should all do fair comparisons, thatâs why I pointed out the issues in your comparisons earlier. I am planning to expand my benchmarking codes to include some more differential equation problems and once I have this new code added with new benchmarking numbers, I will update my results (with and without event callbacks) in my repo. So stay tuned. And you may find some comfort in knowing that for Lorenz system, Julia Tsit5() solver seem to be at least 2x faster than FLINTâs DOP54 (which is Dormand-Prince 5th order or RK45 in Matlab) method in my preliminary testing, but more to come later.
The changes to the callbacks, as I showed above are performance critical. The original code is type unstable and refers to non-const global variables. The new code does not, and as a result is 3x faster.
May I ask perhaps an irrelevant and naive question.
It seems many times if I need to increase Juliaâs speed, I need to use some packages like StaticArray or some macros. Is it possible to just write native Julia code such than I can just reach the speed of using StaticArray?
I guess StaticArray is written in Julia but it uses many macros so that it can make the speed fast for small arrays. However, this optimization process is not easy for regular users, so they make this package as a wrapper.
Wish Juliaâs compiler can automatically do more optimization.
In Fortran usually I do not need to find packages to boost the performance for some basic operations, because the compiler usually did a good job.
I feel the main strength of DifferentialEquations.jl is that it is package of many different algorithms.
I feel if those algorithms are properly written in Fortran, the performance should be similar with Julia.
Many speed difference are coming from the difference in the algorithms. Perhaps, to some extent, a fair comparison might be to use the same algorithm, and compare Julia vs Fortran.
A StaticArray is just Tuple with some arithmetic defined on it for convenience. Arrays are mutable objects and store their length as runtime data, so they will never be 0 overhead. Julia is giving you more safety here, so it takes a little bit of work to make that safety have 0 overhead.
You could relatively easily rewrite the basic solvers in Fortran. That said, Fortran isnât powerful enough to easily deal with things like sparsity of different structures, different input types (eg DoubleFloat, Float16), symbolic simplification, and all the fancy stuff like that. To rival the Julia ecosystem, you would have to impliment a good type system with generics and dispatch, a symbolic algebra library, and basically the rest of a more versatile programming language (Greenspun's tenth rule - Wikipedia).
It was written from the perspective of a Fortran user, dicusses the differences from the user point of view, and at the end shows how to write your own âstatic arraysâ. Fortran indeed does make the user life easier for new users, particularly by optimizing things which, in Julia (still) allocate a new intermediate array at every iteration:
julia> function f!(x)
for i in 1:1000
y = [ x[1], x[2] ] # this allocate a new intermediate array at every iteration
x .= x .+ y
end
return x
end
f! (generic function with 1 method)
julia> @btime f!([1.0,2.0])
27.807 Îźs (1001 allocations: 93.84 KiB)
2-element Vector{Float64}:
1.0715086071862673e301
2.1430172143725346e301
This kind of thing is a frequent cause of very poor first experiences with Julia for people coming from other languages (and one place where a StaticArray, or a tuple, or a custom struct, or just not computing at all the intermediate array, are solutions). That specifically I think that will get optimized at some point in the development of the language (because it is possible to prove that that intermediate array is local to the scope of the loop and, thus, can be stack allocated). But Julia does not do miracles, you have to pave your way into it to write good code in it, as in other languages. I have said this many times, though, in Fortran there is less chances you get a very poor performance to start with. The dynamical character of Julia allows one to write very slow code.
Now, I would suggest you be kind to ask these things on the Julia forum, because it becomes somewhat offtopic here. (also it can be a suggestion for the forum managers here, to split the thread to an âofftopicâ section, I think for future forum users that will be better, particularly because the OP is very important in general for all users).
I have enabled static arrays and some of the changes (Ref values) you are suggesting for event callbacks in my FLINT repo here.
As I had mentioned earlier for this problem (Arenstorf orbit) with or without event callbacks, I am seeing Julia diffeq solvers are slower than FLINT. See my latest results on github.
You (or anyone interested) are welcome to optimize the Julia test script that is used for these tests and submit your changes in my repo directly.
No problem, we all want to see the results more than the claims. I doubt adding const will make it match FLINT numbers as the difference seems to be significant for this problem but I will accept it into the main branch (I think I had removed const originally because it gives those redefinition warnings on repeated runs).
it makes a really big difference because without it being const, globals in Julia are type unstable, which means that you have to do dynamic dispatch, and everything just gets much slower.
You can check the latest benchmarks. I have even wrapped the entire Julia script in a main() function to call from REPL directly instead of including the script as the former method is guaranteed to (JIT) compile everything on the first run.
Thanks for updating the benchmarks! One last question, in the readme it says you are benchmarking vs DifferentialEquations 2.0. If so, you should definitely update to the most recent version (6.19).
Readme probably is not up to date. These latest tests were using DiffEq v6.16. The main thing that improved Julia code performance was when it was run interactively from REPL instead of doing include script every time. Still, the flint solvers can be up to ~5x faster when event detection is on and I think there is still potential for more optimizations in flint code.
One question, is it possible to provide a Makefile,
so perhaps we may just include your test f90 files + all the src f90 files, and do a make, and run?
With the makefile, perhaps we can know the structure/dependence a little easier, once the Makefile works on Linux, I can easier do a makefile for windows, too.
For flint, it is similar, right?
I just need to check your src + test files, and find the dependence and I can write a similar makefile, right? Sorry another lazy question.
Btw, your code looks nice! I can learn a lot from your code! Thank you very much!
Thank you very much!
When you have time, if you could add a makefile it could be great!
It is something like the one in @certik repo,
But anyway, I think I can figure it out. Eg, instead of check the code and do it manually, I can use a IDE software called codeblocks, it can generate a makefile which will show the dependence of each file.
Anyway, great job!
Thank you both @prince_mahajan@oscardssmith
