Puzzle computing arrays

[Tim Prince]

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.