Hi All,
Many ODE solvers use some form of compensated summation to reduce roundoff 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+(tempyi) !line 4
ycsi=ycsi+h*sume
y(i)=yi+ycsi
ycs(i)=ycsi+(yiy(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 longterm errors than the Fortran 90 version (we print out values of conserved quantities, and these numbers drift owing to roundoff 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.108420E18, where 1.0+eps=1.0. My understanding is that somehow they make use of these 80bit 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 roundoff error by roughly a factor of 1000 without a huge decrease in performance.
Jerry