--- Begin Message ---
Tobias Schlüter ha scritto:
Angelo Graziosi wrote:
Tim Prince ha scritto:
Your native code could produce different round-off according to
whether you specified -mfpmath=387 or -mfpmath=sse.
The test case was build with:
gfortran test_case.f90 -o test_case
and the *principle* question still is: why things computed with
implied DO are (a little) different when computed with an explicit DO.
x87 floating point has the annoying property of using 80bit internally,
and only rounding to 64bit when storing values in memory. If we
evaluate at compile-time we always evaluate with REAL*8 or REAL*4
precision. It is hard to predict if the optimizers will end up storing
values to memory, so there's no way in this environment of guaranteeing
if calculations will have the same result if evaluated at compile-time
and if evaluated at runtime or even if two different loops doing the
same calculation will compute the same result. Two ways of making sure
are using SSE-math (which uses 64-bit throughout) or using -ffloat-store
(which stores all values to memory before processing them further).
Now I have understood the nature of the problem. Thanks!
So it would indeed be valuable to see if these options fix your problem.
With -ffloat-store,
$ cat test_case.f90
program test_case
implicit none
integer :: k
integer, parameter :: DP = kind(1.D0),&
N = 3
real(DP), parameter :: A(0:N) = &
(/1.0000000000D0,3.8860009363D0,7.4167881843D0,&
9.8599837415D0/)
real(DP), parameter :: BB(0:N-1) = &
(/((A(k+1)*(A(1)+(k+1))),k = 0,N-2),(A(N)*(A(1)+N))/)
real(DP) :: b(0:N-1),bbb(0:N-1)
do k = 0,N-2
b(k) = A(k+1)*(A(1)+(k+1))
enddo
b(N-1) = A(N)*(A(1)+N)
bbb(0:N-1) = (/((A(k+1)*(A(1)+(k+1))),k = 0,N-2),(A(N)*(A(1)+N))/)
do k = 0,N-1
!print *, k,b(k),BB(k),b(k)-BB(k)
print *, k,b(k)-BB(k),b(k)-bbb(k),BB(k)-bbb(k)
enddo
end program test_case
$ gfortran -ffloat-store test_case.f90 -o test_case
the result is
$ ./test_case.exe
0 0.0000000000000000 0.0000000000000000 0.0000000000000000
1 0.0000000000000000 0.0000000000000000 0.0000000000000000
2 0.0000000000000000 0.0000000000000000 0.0000000000000000
Instead the SSE-math does not work on Cygwin:
$ gfortran -mfpmath=sse test_case.f90 -o test_case
f951: warning: SSE instruction set disabled, using 387 arithmetics
The build of GFortran 4.3.1 I did was
$ gfortran -v
Using built-in specs.
Target: i686-pc-cygwin
Configured with: /work/gcc-4.3.1/configure --prefix=/usr/local/gfortran
--exec-prefix=/usr/local/gfortran --sysconfdir=/usr/local/gfortran/etc
--libdir=/usr/local/gfortran/lib --libexecdir=/usr/local/gfortran/lib
--mandir=/usr/local/gfortran/share/man
--infodir=/usr/local/gfortran/share/info --enable-languages=c,fortran
--enable-bootstrap --enable-decimal-float=bid --enable-libgomp
--enable-threads --enable-sjlj-exceptions
--enable-version-specific-runtime-libs --enable-nls
--enable-checking=release --disable-fixed-point --disable-libmudflap
--disable-shared --disable-win32-registry --with-system-zlib
--without-included-gettext --without-x
Thread model: posix
gcc version 4.3.1 (GCC)
Can I enable SSE? How?
In any case, thanks for clarification.
Cheers,
Angelo.
--- End Message ---