topic After thinking about my post in IntelĀ® Fortran Compiler
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073027#M119789
<P>After thinking about my post #3, it would have worked as intended if N were a power of 2. The following removes any restrictions on N (other than it be positive non-zero and fit in memory).</P>
<PRE class="brush:fortran;">REAL(4) :: T(N/2+1) ! No restrictions on N
REAL(4) :: SUMA
INTEGER(8) :: I,J
...
J = N
I = J / 2
T = A(1:I) + A(I+1:I*2) ! In event J is odd, sum all but last
IF(I*2 .NE J) THEN ! (if N (now J) were odd)
I = I+1
T(I) = T(J) ! copy last to after partial sum
ENDIF
DO WHILE(I .gt. 1)
J = I
I = J / 2
T(1:I/2) = T(1:I/2) + T(I/2+1:I*2)
IF(I*2 .NE J) THEN
I = I+1
T(I) = T(J)
ENDIF
END DO
SUMA = T(1)</PRE>
<P>Performing the same thing on the double arrays would yield a better result.</P>
<P>Jim Dempsey</P>Thu, 19 Jan 2017 13:39:00 GMTjimdempseyatthecove2017-01-19T13:39:00ZBug in intrinsic function sum() compiled with -fp-model precise
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073023#M119785
<P>I'm running the latest ifort 2017.1 and found this unwanted behavior. It's illustrated by a short program and the result after two different compilations.</P>
<P>> cat test_sum.f<BR />
PROGRAM TEST_SUM<BR />
INTEGER(8),PARAMETER :: N = 1024 * 1024 * 1024<BR />
INTEGER(8) :: SEED(1) = (/12345678/)<BR />
REAL(4) :: A(N)<BR />
REAL(8) :: B(N)</P>
<P>! Fill the array A with random numbers in the range [0 - 1].<BR />
CALL RANDOM_SEED(PUT=SEED)<BR />
CALL RANDOM_NUMBER( A )</P>
<P>! Copy the values to the double precision array B.<BR />
B(:) = A(:)</P>
<P>! Expected sum = N * 0.5 = 536870912 = 5.3687091E+08<BR />
! The sum may differ a bit due to the random numbers drawn.<BR />
PRINT *, 'Sum(A):', REAL(SUM(A))<BR />
PRINT *, 'Sum(B):', REAL(SUM(B))</P>
<P> END PROGRAM TEST_SUM<BR />
><BR />
> ifort -i8 test_sum.f<BR />
> a.out<BR />
Sum(A): 5.3684355E+08<BR />
Sum(B): 5.3687296E+08<BR />
><BR />
> ifort -i8 -fp-model precise test_sum.f<BR />
> a.out<BR />
Sum(A): 1.6777216E+07<BR />
Sum(B): 5.3687296E+08</P>
<P>So the real*8 array was summed up correctly, but not the real*4 array when compiled with the -fp-model precise option. The result was the same if the arrays A and B were allocated arrays.</P>Wed, 18 Jan 2017 05:54:19 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073023#M119785Lars-goran_Ofversted2017-01-18T05:54:19Z>>Expected sum = N * 0.5 =
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073024#M119786
<P>>>Expected sum = N * 0.5 = 536870912 = 5.3687091E+08</P>
<P>Worst case difference could be N lsb's as you accumulate each partial sum. Average may be N*0.5 lsb's</P>
<P>Float/Single has 7.22 digits of precision. 10^7.22 = 16595869 ~= 1.66E+7</P>
<P>When the partial sum of the float approaches and exceeds 1.66E7, the additional additions effectively add 0.0.<BR />
The summation of the float array is off by (5.3687091E+08 - 1.66E+7) ~= 5.2E+8. The actual difference will depend on when the partial sum approaches and exceeds 8.3E+6. The aove is thumbnail sketch, a mathematician can provide a proper statistical analysis.</P>
<P>The results appear realistic.</P>
<P>Jim Dempsey</P>Wed, 18 Jan 2017 18:03:54 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073024#M119786jimdempseyatthecove2017-01-18T18:03:54ZTry merging this into your
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073025#M119787
<P>Try merging this into your program and report back the results:</P>
<PRE class="brush:fortran;">REAL(4) :: T(N/2) ! N must be even (it is)
REAL(4) :: SUMA
INTEGER(8) :: I
...
T = A(1:N/2) + A(N/2+1:N)
I = N / 2
DO WHILE(I .gt. 1)
T(1:I/2) = T(1:I/2) + T(I/2+1:I)
I = I / 2
END DO
SUMA = T(1)</PRE>
<P>Jim Dempsey</P>Wed, 18 Jan 2017 18:17:05 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073025#M119787jimdempseyatthecove2017-01-18T18:17:05ZIn case someone wonders why
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073026#M119788
<P>In case someone wonders why -fp-precise might be responsible for such a difference, that option requires SUM to be implemented as a single chain of serial addition. Without that option, ifort can choose to optimize performance [ by splitting the sum into several partial sums (typically 8 individual 4- or 8- wide parallel sums, depending on target ISA), and this would delay the effect Jim described.</P>Wed, 18 Jan 2017 19:39:22 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073026#M119788TimP2017-01-18T19:39:22ZAfter thinking about my post
https://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073027#M119789
<P>After thinking about my post #3, it would have worked as intended if N were a power of 2. The following removes any restrictions on N (other than it be positive non-zero and fit in memory).</P>
<PRE class="brush:fortran;">REAL(4) :: T(N/2+1) ! No restrictions on N
REAL(4) :: SUMA
INTEGER(8) :: I,J
...
J = N
I = J / 2
T = A(1:I) + A(I+1:I*2) ! In event J is odd, sum all but last
IF(I*2 .NE J) THEN ! (if N (now J) were odd)
I = I+1
T(I) = T(J) ! copy last to after partial sum
ENDIF
DO WHILE(I .gt. 1)
J = I
I = J / 2
T(1:I/2) = T(1:I/2) + T(I/2+1:I*2)
IF(I*2 .NE J) THEN
I = I+1
T(I) = T(J)
ENDIF
END DO
SUMA = T(1)</PRE>
<P>Performing the same thing on the double arrays would yield a better result.</P>
<P>Jim Dempsey</P>Thu, 19 Jan 2017 13:39:00 GMThttps://community.intel.com/t5/Intel-Fortran-Compiler/Bug-in-intrinsic-function-sum-compiled-with-fp-model-precise/m-p/1073027#M119789jimdempseyatthecove2017-01-19T13:39:00Z