Poisson 2D Jacobi solver benchmark with Python, Cython, C, and Fortran code

[arunningcroc]

As discussed in the discourse group.

[certik]

Thank you @arunningcroc!

@milancurcic, @LKedward, @awvwgk, how should we proceed?

My goal is to just start adding benchmarks, even if we don't (yet) have an infrastructure for it, and then simply later adding some tools to run them easily (with different compilers and options), etc.

[smeskos]

Hey all, long time no see... If I may, I have three remarks:

1. At Readme.md change the sentence: And the timing results (the grid is MxM): to something like "The grid is a matrix of size MxM and the timing results in seconds: and then remove (seconds) from within the table's header row. It was confusing to me. Also, add the symbol * before the sentence below the table.
2. Perhaps gfortran optimizes it automatically but I would change the loop order to reflect the Column-major aspect of it, even just for the appearances.
3. There is an inadvertent error at the optimized.c -> swap function, you just assigned B to A, instead of swapping.

[milancurcic]

Thank you @arunningcroc!

My suggestion is to add only Fortran code, until we have a clear plan on whether, why, and how to include other languages.

[certik]

We discussed quite a few times, e.g. in #2, to include other languages.

Why: to actually have a comparison across languages. I expect that across Fortran/C/C++ perhaps even Julia one can fiddle with the code to eventually get similar speed. And we should have that final code. But also what would be interesting (for me) would be code that you would reasonably write as a domain expert, say a physicist. And we should have those codes too.

How: that has to be discussed, for now let's have it in some form, such as in this PR. I expect we'll have benchmarks for smaller arrays, that have to be run many times. I have a runner for that somewhere, I'll see if I can contribute it. Then I expect to have longer running tests, such as this PR, which are not as sensitive to how it is timed.

[milancurcic]

I opened #10.

[milancurcic]

I left a few suggestions in formatting the table.

We should also note what hardware and OS was this run on.

I suggest .f95 -> .f90, at least until we reclaim .f.

[certik]

This looks good to me to merge.

I am going to wait for @milancurcic to be ok with it.

[milancurcic]

I approved earlier, I think it's great. Thank you!

[certik]

@arunningcroc, we go some feedback on the code. When compiling with gfortran -Ofast optimized.f90 and decreasing the problem size to M=100 yields non-deterministic answers. Is this expected? Is there some use of uninitialized data?

There should be a verification phase (not part of the benchmark) to ensure the result is actually correct, and that will help when people writing alternative implementations to check that their code is correct. Preferably a solution that works for all M.

The final feedback we got was that the documentation could be improved in README, it would be good to document the numerical method a bit, the boundary/initial condition, etc.

@arunningcroc I can help with this. Let me know.

[rouson]

Should this pull request be reopened? I can't do it.

I prompted @certik's most recent comment so I'm pasting a transcript below showing the non-deterministic output. I didn't use any optimization in the transcript below, but -Ofast yields similar results.

Regarding documentation, I suggest one short comment per "major" construct, wherein my personal definition of "major" includes modules; procedures or procedure interfaces, and derived types.

± git diff
diff --git a/poisson2d/optimized.f90 b/poisson2d/optimized.f90
index 3bc1045..316ba64 100644
--- a/poisson2d/optimized.f90
+++ b/poisson2d/optimized.f90
@@ -19,7 +19,7 @@ end module
 program poisson
 use rhofunc, only: rho
 implicit none
-integer, parameter :: dp=kind(0.d0), M=300
+integer, parameter :: dp=kind(0.d0), M=100
 integer            :: i,j, iter
 real(dp),parameter :: epsilon0=8.85E-12_dp, target=1E-6_dp, a=0.01
 real(dp)           :: delta, b, e, phiprime(M,M), phi(M,M), a2, rhoarr(M,M)

± gfortran optimized.f90 
± ./a.out
   1.6580310000000000            22050
± ./a.out
   1.6223630000000000            22050
± ./a.out
   1.7145189999999999            22050

[rouson]

Never mind the issue of non-deterministic results! I just realized the that first number is simply the CPU time. I really hope that submitted codes will use more obvious variable names and add even just a small number of comments in addition to providing results-verification.

[certik]

We could revert master and open a new PR, but I wouldn't bother. But we'll fix this, even if I have to do it myself. But if @arunningcroc has time, he can help me.