Puzzle computing arrays

[Angelo Graziosi]

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.
> 

So, which is the best solution?

Your conclusion is not clear to me, sorry!


Cheers,
    Angelo.