Ifort (IFORT) 2021.8 evaluates `1.0E+37 / 1.0E+38` to 0

[Zaikun]

This is related to the interesting behavior of ifort discussed before, but the result is even stranger. 

N.B.: Only the lastest ifort, namely Ifort (IFORT) 2021.8.0, behaves in this way. ifort 2021.7.1 or ifx does not.

The same question is posed on Fortran Discourse.

Code:

program test

implicit none

real :: S(7, 7)

S = 0.0
S(1, 7) = 1.0E+37
S(7, 1) = S(1, 7)

write (*, *) 'Before: S(1, 7), S(7, 1)', S(1, 7), S(7, 1)

! Is S symmetric?
if (all(abs(S - transpose(S)) <= 0)) then
    write (*, *) 'S is symmetric.'
else
    error stop 'S is not symmetric, which is wrong'
end if

S = S / 1.0E+38

!write (*, *) S
write (*, *) 'After: S(1, 7), S(7, 1)', S(1, 7), S(7, 1)

! Is S symmetric?
if (all(abs(S - transpose(S)) <= 0)) then
    write (*, *) 'S is symmetric.'
else
    error stop 'S is not symmetric, which is wrong'
end if

end program test

Result:

$ ifort --version && ifort test_intel_sym.f90 && ./a.out  
ifort (IFORT) 2021.8.0 20221119
Copyright (C) 1985-2022 Intel Corporation.  All rights reserved.

 Before: S(1, 7), S(7, 1)  9.9999999E+36  9.9999999E+36
 S is symmetric.
 After: S(1, 7), S(7, 1)  0.0000000E+00  0.1000000 
 S is not symmetric, which is wrong

What is happening? Thanks.

[Zaikun]

Update: I suppose this is the simplest version of the test.

program test
implicit none

integer :: i
real :: x(5)
x = 1.0E37  ! All elements of x are set to the same value.
print *, 'Before division, x = ', x
print *, 'All entries of x are identical?', all(x == x(1))
x = x / 1.0E38 
print *, 'After division, x = ', x
do i = 1, size(x)
    print *, 'i = ', i, 'x(1) == x(i)?', x(1) == x(i)

    print "(*(b0,1x))", x(1), x(i)

    if (x(1) /= x(i)) then
        print *, ''
        print *, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        print *, '!!!!!!!!!!! Surprise! x(1) is different from x(i) !!!!!!!!!'
        print *, '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
        print *, ''
    end if
end do

end program test

$ uname -a && ifort --version && ifort test_div.f90  && ./a.out
Linux zX10 5.15.0-56-generic #62-Ubuntu SMP Tue Nov 22 19:54:14 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux
ifort (IFORT) 2021.8.0 20221119
Copyright (C) 1985-2022 Intel Corporation.  All rights reserved.

 Before division, x =   9.9999999E+36  9.9999999E+36  9.9999999E+36
  9.9999999E+36  9.9999999E+36
 All entries of x are identical? T
 After division, x =   0.0000000E+00  0.0000000E+00  0.0000000E+00   0.0000000E+00  0.1000000    
 i =            1 x(1) == x(i)? T
0 0
 i =            2 x(1) == x(i)? T
0 0
 i =            3 x(1) == x(i)? T
0 0
 i =            4 x(1) == x(i)? T
0 0
 i =            5 x(1) == x(i)? F
0 111101110011001100110011001101
 
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 !!!!!!!!!!! Surprise! x(1) is different from x(i) !!!!!!!!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

[Umar__Sait]

It seems like if you declare x and the exponential constants as double precision it works.

[Zaikun]

Right, but that does not resolve the problem. Thank you. 

[FortranFan]

@Zaikun ,

You may want to take note your posts are getting repetitive and unrevealing of any new issues.

You were informed weeks ago to apply -fp:strict with the floating-point computational instructions on interest to you with Intel Fortran, especially if you are working with the default REAL intrinsic type in your Fortran code.

I doubt you will be able to reproduce the problem here listed in your original post here or the second one thereafter with -fp:strict.  If you do, then post back - otherwise note what you are bringing up are simply how things are, basically non-issues, and they are so for the reasons pointed out to you earlier.

[Zaikun]

Thank you for pointing things out. I am sorry that I am not posting something that interests you.

I can understand that `S(1,7)` differs from `S(7,1)` due to vectorization and lack of `-fp-model=strict`.

However, I do not believe that fp model is a sufficiently good excuse for evaluating `1.0E37 / 1.0E38` to `0`.

If it is quite desirable to have `1.0E37 / 1.0E38 = 0`, why don't `ifx` and earlier versions of `ifort` behave in this way?

I am afraid that it is a bug.

Thanks.

[Ron_Green]

I will investigate this one @Zaikun 

[FortranFan]

@Ron_Green ,

Thank you for offering to OP to investigate this.

Please note this is a case where OP and other such users are getting "astonished", or suffering from worse reactions.

Please consider Intel's documentation for `-prec-div`:

https://www.intel.com/content/www/us/en/develop/documentation/fortran-compiler-oneapi-dev-guide-and-reference/top/compiler-reference/compiler-options/floating-point-options/prec-div-qprec-div.html

It says clearly, "At default optimization levels, the compiler may change floating-point division computations into multiplication by the reciprocal of the denominator. For example, A/B is computed as A * (1/B) to improve the speed of the computation."

Please note however there is no mention of vectorization around the compiler optimization nor about anything regarding the precision of the floating-point type.  As such, a user might expect similar issues with scalar division or those involving whole arrays and for the approximate result to be present whether the floating-point type is of single, or double precision, etc.  However please also consider the case below:

   integer, parameter :: WP = selected_real_kind( p=6 )
   real(WP), parameter :: BIG = 1.0_wp / tiny( 1.0_wp )
   real(WP), parameter :: FOO = 0.12_wp*BIG
   real(WP), parameter :: BAR = 1.2_wp*BIG
   real(WP) :: x(3)
   blk1: block
      print *, "Block 1: unrolled division"
      x = FOO
      print *, "x = ", x
      x(1) = x(1) / BAR ; x(2) = x(2) / BAR ; x(3) = x(3) / BAR
      print *, "After division: x = ", x
   end block blk1
   print *
   blk2: block
      print *, "Block 2: vector division"
      x = FOO
      print *, "x = ", x
      x = x / BAR 
      print *, "After division: x = ", x
   end block blk2
end

C:\temp>ifort /standard-semantics p.f90
Intel(R) Fortran Intel(R) 64 Compiler Classic for applications running on Intel(R) 64, Version 2021.7.0 Build 20220726_000000
Copyright (C) 1985-2022 Intel Corporation.  All rights reserved.

Microsoft (R) Incremental Linker Version 14.33.31630.0
Copyright (C) Microsoft Corporation.  All rights reserved.

-out:p.exe
-subsystem:console
p.obj

C:\temp>p.exe
 Block 1: unrolled division
 x =  1.0208471E+37 1.0208471E+37 1.0208471E+37
 After division: x =  0.1000000 0.1000000 0.1000000

 Block 2: vector division
 x =  1.0208471E+37 1.0208471E+37 1.0208471E+37
 After division: x =  0.000000 0.000000 0.1000000

C:\temp>

Clearly several artifacts of Intel's optimization around vectorization play a role with the case in "Block 2" where ostensibly the first two elements of the array 'x' are affected but not the third.

From a product stewardship point-of-view with Intel Fortran, the program behavior in the case shown here can be accepted as a violation of the "principle of least astonishment" when it comes to software design.  It will be good if you can discuss this with your team and provide comments.

[jimdempseyatthecove]

>>lack of `-fp-model=strict`....However, I do not believe that fp model is a sufficiently good excuse for evaluating `1.0E37 / 1.0E38` to `0`.

Create a configuration file containing -fp-model=strict (and any additional options you see fit for default behavior).

The problem is when (the default) -fp-model=fast is used, some divisions are performed via the faster two-step process of generate an approximate reciprocal then multiply. You need to instruct the compiler to "don't do that".

Jim Dempsey

[jimdempseyatthecove]

Perhaps Intel should consider a firmware change to the reciprocal instruction whereby whenever the exponent nears its limit (E+../E-..) that formal (long) division is performed to produce the more accurate reciprocal value.

Jim Dempsey

[Umar__Sait]

I agree that something needs to be done. I teach Fortran to college and graduate students who would be very confused with this behavior, specially astronomy students who use everything in cgs units with large exponents. Somehow the under the hood optimization tricks should be more transparent to the user.

[Ron_Green]

I opened bug ID CMPLRIL0-35107 for the recip-mult and vectorization concerns. I will open a DOCumentation issue to add more warnings to prec-div

[Zaikun]

Thank you @Ron_Green for taking care of this. 

[Ron_Green]

@FortranFan the documentation task ID is DOC-10684

[FortranFan]

Thank you @Ron_Green , that will be very helpful for users generally.