This is the mail archive of the fortran@gcc.gnu.org mailing list for the GNU Fortran project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Re: Puzzle computing arrays


Angelo Graziosi wrote:
Tim Prince ha scritto:
Angelo Graziosi wrote:
The following happens on Cygwin with GFortran 4.3.1 and on GNU/Linux Kubuntu 8.04 with GFortran 4.3.1 and 4.2.3.

In the following test case (reduced to the simplest form), the array elements of BB() and b() are computed exactly with the same expression (copy/pasted), but the results seem a little different (only on 14-15th digit)

$ cat test_case.f90
!
! gfortran test_case.f90 -o test_case
!
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)
  do k = 0,N-2
     b(k) = A(k+1)*(A(1)+k+1)
  enddo
  b(N-1) = 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)
  enddo
end program test_case

$ gfortran test_case.f90 -o test_case

$ ./test_case.exe
           0   0.0000000000000000
           1 -7.10542735760100186E-015
           2   0.0000000000000000

Since BB() and b() are computed with the same expression, I would expect the same results. But BB() is a PARAMETER (computed at compile time???) and b() is computed at run time...

So, what could be the difference? What explanation may you give?

Don't the compile time calculations get extra precision from mpfr and friends? The run-time calculation would make more sense with parentheses (k+1) to specify a single round-off.


Another test case (the previous++):


$ cat test_case.f90
!
! gfortran test_case.f90 -o test_case
!
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

$ ./test_case.exe
0  0.0000000000000000        0.0000000000000000       0.0000000000000000
1 -7.10542735760100186E-015 -7.10542735760100186E-015 0.0000000000000000
2  0.0000000000000000        0.0000000000000000       0.0000000000000000

Now bbb() is NOT a PARAMETER but is still computed with an implied DO.

It looks that computing the same expression with an implied DO is a little different when using an explicit DO.

Ideally, a short array constructor would be expanded at compile time, when that is possible, using mpfr rather than native code. Your native code could produce different round-off according to whether you specified -mfpmath=387 or -mfpmath=sse.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]