Many ODE solvers use some form of compensated summation to reduce round-off errors. In my particular case, there is a code segment that looks like this:
temp=y(i) ycsi=ycs(i)+h*suma yi=temp+ycsi ycsi=ycsi+(temp-yi) !line 4 ycsi=ycsi+h*sume y(i)=yi+ycsi ycs(i)=ycsi+(yi-y(i)) !line 7
The order of operations is important for compensated summation to work. For example, one might manipulate the above expressions and conclude that ysci at line 4 is zero, and that ysc(i) at line 7 is zero.
We have this ODE solver implemented in Fortran 90 and also Matlab. The Matlab version has slightly smaller long-term errors than the Fortran 90 version (we print out values of conserved quantities, and these numbers drift owing to round-off error). I strongly suspect that the Fortran compiler is tinkering with the code, where Matlab does not.
How can I tell the compiler not to tinker with that bit of code? I tried using -O0 at compile time, but that did not seem to make any difference.
Better yet, is there a compiler independent way to do this? For example, one can use OpenMP directives and many (if not most) Fortran compilers understand these.
I noticed that gfortran has a REAL*10 (KIND=10) option. The machine precision for this option is eps = 0.108420E-18, where 1.0+eps=1.0. My understanding is that somehow they make use of these 80-bit thingamabobs that most Intel CPUs have, and as a result, the speed difference is only about a factor of 2 slower than when using KIND=8. Are there plans to include KIND=10 in the PGI compilers? In the case of our ODE solver, going to KIND=10 beats down the round-off error by roughly a factor of 1000 without a huge decrease in performance.