------- Comment #3 from dominiq at lps dot ens dot fr 2007-10-15 15:51 -------
> that causes numerical results with CP2K to change going from -O0 to -O3.
If you do expect that optimization optimizes your computation, you should
expect some change of the numerical results, so put some tolerance in the
comparisons.
How small can be the tolerance depends on your problem. I have modified your
code to compute the polynomial in three different ways (yours, plus two others,
the third one being the recommended one for "numerical stability", if I did not
make some errors):
FUNCTION F(r,a,e,pf,qf)
REAL*8 :: a(2:15)
REAL*8 :: r,f,e, pr,pf,qf,rr
f=0
pf=0
pr=r
rr=1.0d0/r
DO i = 2, 15
f = f + a(i)/(r**(i-1)*REAL(i-1,8))
pf = pf + a(i)/(pr*REAL(i-1,8))
pr = r*pr
END DO
f=f-e
pf=pf-e
qf=a(15)/REAL(14,8)
do i = 14, 2, -1
qf=rr*qf+a(i)/REAL(i-1,8)
end do
qf=rr*qf-e
END FUNCTION F
PROGRAM TEST
REAL*8 :: a(2:15)=(/-195.771601327700D0, 15343.7861339500D0, &
-530864.458651600D0, 10707934.3905800D0, &
-140099704.789000D0, 1250943273.78500D0, &
-7795458330.67600D0, 33955897217.3100D0, &
-101135640744.000D0, 193107995718.700D0, &
-193440560940.000D0, -4224406093.91800D0, &
217192386506.500D0, -157581228915.500D0/)
REAL*8 :: r=4.51556282882533D0
REAL*8 :: e=-2.21199635966809471000D0
REAL*8 :: pf, qf
REAL*8 :: f
write(6,*) f(r,a,e,pf,qf)
write(6,*) pf, qf
END PROGRAM TEST
[karma] f90/bug% gfc pr33780_red.f90
[karma] f90/bug% a.out
0.0000000000000000
1.36957112317759311E-012 7.66053886991358013E-012
[karma] f90/bug% gfc -O3 pr33780_red.f90
[karma] f90/bug% a.out
-1.14397380457376130E-012
1.36957112317759311E-012 1.04861186749364478E-011
So for your problem the tolerance should be ~1.0E-11, and since pf and qf don't
use powers, this is not a question of inlining, but how the operations are
shuffled to "optimize" your computation.
--
http://gcc.gnu.org/bugzilla/show_bug.cgi?id=33780
