Actually, I ran it again a few times, and now Fortran seems faster. Now it seems to converge in 320k iterations (for 250x250) compared to C’s 304k, but is still faster by a bit less than a second. Also, I used epsilon = 8.85E-12, with the smaller epsilon it converges in 100k less iterations and is then obviously much faster.
Also, the C code becomes a few seconds slower if you allocate on the heap. Not sure how Fortran handles that behind the scenes (big arrays without allocate), but for a realistic numerical program you’d probably want to malloc in C.
I’d be willing to give it a shot, if there are some guidelines on how to write these benchmarks. I also have a Cython version of this code, as well as pure python, and it would of course trivial to write a NumPy array slicing implementation.
Also, I got the looping Fortran version to run in 15 seconds now for 250x250, by following @implicitall 's lead and removing the useless temp array as well as precalculating rho array.This is now some 4 seconds faster than C with malloc.
We don’t have those yet, but I think it’s better to simply get started and go from there. Why don’t you submit a PR against the benchmarks repository, just create some directory and put the benchmarks there. Thanks!
Just a word of caution. Using -Ofast can have unexpected compute results in certain cases (loss of accuracy). I never use it. I generally use -O2 and sometimes -O3. Speed is not everything.