Strange behavior of `ifort`

1

Consider the following code.
[Update: @septc has simplified the code to avoid the strange numbers]

! test.f90                                                                                                                           
program test                                                                                                                         
                                                                                                                                     
implicit none                                                                                                                                                                                                        
                                                                                                                                     
! Why these strange numbers? Because they come from a real project.                                                                  
real, parameter :: A(25) = &                                                                                                         
    & [1.2409463E+22, 2.8192031E+20, -4.4971432E+21, 1.1558917E+22, -4.7725394E+21, &                                                
    & 2.8192031E+20, 2.8192031E+20, -3.2276437E+21, -3.3936284E+21, 4.1412172E+21, &                                                 
    & -4.4971432E+21, -3.2276437E+21, -7.6997325E+20, -4.1627641E+21, -3.4729614E+20,&                                               
    & 1.1558917E+22, -3.3936284E+21, -4.1627641E+21, 2.5757004E+21, -2.4898961E+21, &                                                
    & -4.7725394E+21, 4.1412172E+21, -3.4729614E+20, -2.4898961E+21, 3.5042103E+21]                                                  
                                                                                                                                     
real :: S(5, 5)                                                                                                                      
                                                                                                                                     
S = reshape(A, [5, 5])                                                                                                               
                                                                                                                                     
! Is S symmetric?                                                                                                                    
write (*, *) all(abs(S - transpose(S)) <= 0)                                                                                         
                                                                                                                                     
! The numbers in S are horrible. Let's scale it.                                                                                     
S = S / maxval(abs(S))                                                                                                               
                                                                                                                                     
write (*, *) S                                                                                                                       
                                                                                                                                     
! Is S symmetric?                                                                                                                    
write (*, *) all(abs(S - transpose(S)) <= 0)                                                                                         
                                                                                                                                     
end program test

Here is what’s happening on my computer.

$ uname -a && ifort --version 

Linux 5.15.0-53-generic #59-Ubuntu SMP Mon Oct 17 18:53:30 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

ifort (IFORT) 2021.7.1 20221019
Copyright (C) 1985-2022 Intel Corporation.  All rights reserved.

$ ifort test.f90 && ./a.out 
 
T

  0.9999999      2.2718171E-02 -0.3623963      0.9314599     -0.3845887    
  2.2718171E-02  2.2718171E-02 -0.2600954     -0.2734710      0.3337145    
 -0.3623963     -0.2600954     -6.2047265E-02 -0.3354508     -2.7986396E-02
  0.9314599     -0.2734710     -0.3354508      0.2075594     -0.2006450    
 -0.3845887      0.3337145     -2.7986394E-02 -0.2006450      0.2823821    

 F

Note that S(3,5) /= S(5,3) after the scaling. This happens when ifort is invoked with -O, -O2, -O3, -O4, or -O5, but not -g, -O0, -O1, or -Ofast.

Is this expected?

3 Likes
2

Not expected, but not unexpected. You are using the default real kind, which is a 32-bit value.
You have specified more digits than are representable. So you can only expect the precision to be within that range (see the PRECISION() and DIGITS() intrinsics). So you basically have to consider digits past the precision of the computation to be unreliable.

So it is a little surprising, but the answers are still the same within the precision you chose.
I could not reproduce it on my platform using ifort, which is interesting. Could you print S
before you manipulate it? You might be surprised the values not exactly representable will not exactly match what you entered. Is that unexpected for you?

3

Hello @urbanjost , thank you for the response. Here is the revised code and the result, with S printed.

! test.f90
program test

implicit none

! Why these strange numbers? Because they come from a real project.
real, parameter :: A(25) = &
    & [1.2409463E+22, 2.8192031E+20, -4.4971432E+21, 1.1558917E+22, -4.7725394E+21, &
    & 2.8192031E+20, 2.8192031E+20, -3.2276437E+21, -3.3936284E+21, 4.1412172E+21, &
    & -4.4971432E+21, -3.2276437E+21, -7.6997325E+20, -4.1627641E+21, -3.4729614E+20,&
    & 1.1558917E+22, -3.3936284E+21, -4.1627641E+21, 2.5757004E+21, -2.4898961E+21, &
    & -4.7725394E+21, 4.1412172E+21, -3.4729614E+20, -2.4898961E+21, 3.5042103E+21]

real :: S(5, 5)

S = reshape(A, [5, 5])

write (*, *) 'S before scaling:'
write (*, *) S

write (*, *) 'S - transpose(S) before scaling:'
write (*, *) S - transpose(S)

! Is S symmetric?
write (*, *) 'Is S symmetric?', all(S == transpose(S))
write (*, *) 'Does S - transpose(S) contain only zero?', all(abs(S - transpose(S)) <= 0)

! The numbers in S are horrible. Let's scale it.
S = S / maxval(abs(S))

write (*, *) 'S after scaling:'
write (*, *) S

write (*, *) 'S - transpose(S) after scaling:'
write (*, *) S - transpose(S)

! Is S symmetric?
write (*, *) 'Is S symmetric?', all(S == transpose(S))
write (*, *) 'Does S - transpose(S) contain only zero?', all(abs(S - transpose(S)) <= 0)

end program test

$ uname -a && ifort --version 

Linux 5.15.0-53-generic #59-Ubuntu SMP Mon Oct 17 18:53:30 UTC 2022 x86_64 x86_64 x86_64 GNU/Linux

ifort (IFORT) 2021.7.1 20221019
Copyright (C) 1985-2022 Intel Corporation.  All rights reserved.

$ ifort test.f90 && ./a.out 
 
 S before scaling:
  1.2409463E+22  2.8192031E+20 -4.4971432E+21  1.1558917E+22 -4.7725394E+21
  2.8192031E+20  2.8192031E+20 -3.2276437E+21 -3.3936284E+21  4.1412172E+21
 -4.4971432E+21 -3.2276437E+21 -7.6997325E+20 -4.1627641E+21 -3.4729614E+20
  1.1558917E+22 -3.3936284E+21 -4.1627641E+21  2.5757004E+21 -2.4898961E+21
 -4.7725394E+21  4.1412172E+21 -3.4729614E+20 -2.4898961E+21  3.5042103E+21
 
S - transpose(S) before scaling:
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00

 Is S symmetric? T
 Does S - transpose(S) contain only zero? T

 S after scaling:
  0.9999999      2.2718171E-02 -0.3623963      0.9314599     -0.3845887    
  2.2718171E-02  2.2718171E-02 -0.2600954     -0.2734710      0.3337145    
 -0.3623963     -0.2600954     -6.2047265E-02 -0.3354508     -2.7986396E-02
  0.9314599     -0.2734710     -0.3354508      0.2075594     -0.2006450    
 -0.3845887      0.3337145     -2.7986394E-02 -0.2006450      0.2823821    

 S - transpose(S) after scaling:
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00 -1.8626451E-09
  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00  1.8626451E-09  0.0000000E+00  0.0000000E+00

 Is S symmetric? F
 Does S - transpose(S) contain only zero? F
4

How about printing S- transpose(S)?

5

  1. What you said is totally understandable to me. However, according to the test all(abs(S - transpose(S)) <= 0) before the scaling, S is indeed symmetric, with the values represented in the memory.

  2. As you can probably conceive, I do not really enter these values in the real project. Instead, the real-project code looks like

S = (some code that is supposed to generate a symmetric matrix) 
                                                                                                                                     
! Is S symmetric?                                                                                                                    
write (*, *) all(abs(S - transpose(S)) <= 0)   ! The result is `TRUE`                                                                                 
                                                                                                                                     
! The numbers in S are horrible. Let's scale it.                                                                                     
S = S / maxval(abs(S))                                                                                                               
                                                                                                                                     
write (*, *) S                                                                                                                       
                                                                                                                                     
! Is S symmetric?                                                                                                                    
write (*, *) all(abs(S - transpose(S)) <= 0)   ! The result is `FALSE`
6

Good idea! To reduce the amount of code in this thread, I updated the response to @urbanjost just now and included it. Thank you for checking.

7

I ran it myself on Windows: with default options ifort gives a difference of 1.8e-9, really in the noise range - 10 million times smaller than the original value of -2.8e-2. Curilously, the two non-zero values are equal in magnitude but opposite signs.

If I use the option -Od, the result is exactly zero.

It might be better to use epsilon() instead of 0 as the bound.

8

Hello,
This is a somewhat simplified version of the code, where I just used 1/real(i) to generate array elements and made the array S symmetric afterward. Then, if I divide S by tmp, the S becomes not symmetric. If I attach the -fp-model stric or -check all option (in addition to -O3), the S becomes symmetric again. The result also becomes symmetric if I multiply by tmp (rather than division).

program test                                                                            
implicit none                                                                                                                                                                                                        
integer :: i                                                                                                                              
real, parameter :: A(25) = [( 1 / real(i), i=1,25 )]
real :: S(5, 5), tmp                                                                                                           
                                                                                                                                     
S = reshape(A, [5, 5])                                                                                                               
S = ( S + transpose(S) ) / 2

print *, "S(before) = "
print *, S
print *, "symmetric? ", all(abs(S - transpose(S)) <= 0)                                                                                         

tmp = 10                                                                                                                                     
S = S / tmp   !! not symmetric (-O3), symmetric (-O3 -fp-model strict or -chek all)
!! S = S * tmp     !! symmetric (-O3)

print *
print *, "S(after) = "
print *, S                                                                                                                                    
print *, "symmetric? ", all(abs(S - transpose(S)) <= 0)      

print *
print *, "error in symmetry"
print *, abs(S - transpose(S))
                                                                                                                                     
end program

I am not familiar with recent ifort options, but apart from the precision stuff, is it possibly related to some optimization that breaks the symmetry of S…?

(I tried it with Compiler explorer here. code + result )

2 Likes
9

Does the program claim that S is symmetric before the scaling? If yes, is it expected/acceptable that S becomes asymmetric after a floating point division by a scalar? IMHO, the “noise” in floating point arithmetic should not behave in such a way.

This is not surprising. Similarly, the behavior mentioned in my post manifests itself only with ‘-O2’ or above, but not with ‘-g’, ‘-Ofast’, ‘-O0’, or ‘-O1’ on Linux.

10

Thank you @septc very much for this simplified code!

To be totally frank, I do not think any reasonable “optimization” should lead to such behavior. An optimization should at least have some benefit, yet it is difficult to imagine what would the benefit of doing so. It is a bug, or at least a bizarre “feature”.

11

Yeah, I initially imagined that it is related to just floating-point things, but after trying the above code (with tmp), I also felt something strange (because it is just a uniform scaling of all elements). For comparison, I also tried ifx and gfortran (both the newest versions) with -O3 or -Ofast, but they give symmetric S. I am not sure this is also something not too surprising as a result of optimization (including floating-point things)…

1 Like
12

On MacOS, ifort 2021.7.1 all non-zero values of abs(S-transpose(S)) in @septc version are 4.6566129E-10 but this is not a “common” value. When output in hex, it shows to be just a single bit: 0X8.00000000P-34

1 Like
13

I can confirm this issue for ifort.

However, if I define

integer, parameter :: r8=selected_real_kind(15,9)

and use

real(kind=r8)

instead of real. Then it print T and T.

Perhaps can go to Intel website and report this issue.

14

Did you try

write (*, *) all(abs(S - transpose(S)) <= epsilon(1.0) )

with the original? Something like that, as suggested above, is a superior test in general; there are some general guides listed on the fortran.lang and Fortran Wiki site that are a good read about floating point precision. Even with an epsilon(1.0) check if you still get an error that will be significant to report to the Intel forum.

1 Like
15

+1 :+1:

16

I think the issue here is that the quantity x/y is computed twice, within the same array expression in fact, and two different values are produced. Even with inaccurate floating point arithmetic, that should not happen, one should at least expect reproducible results. The difference is in the least significant bit, so it looks like a rounding mode flag is being reset or corrupted somehow within that array expression.

1 Like
17

After eliminating the array syntax and using DO loops it got even odder, as these two
loops did not produce the same answer; where the second one had consistent results;

S2=S
big=maxval(abs(S))

do j=1,5,1 
do i=1,5,1 
      S(i,j)=S(i,j)/big
   enddo
enddo

do j=5,1,-1
do i=5,1,-1
      S2(i,j)=S2(i,j)/big
   enddo
enddo

and at the top I printed everything with a b32.32 format and the elements were identical, but
before computation if I did “S=nearest(S)” then the results became consistent in the original.
Definitely seems like Intel should explain this one, as even if looked as *.s files and so on cannot change it anyway. Have some wild guesses on what is happening, so will be interested in what is really causing it, though.

1 Like
18

Glad to see that the discussion is on the right track finally.

To be brief, the problem is not about how real numbers are represented and calculated on computers, or how we should compare these numbers. All of these topics are important, but unfortunately, miss the point under the current discussion.

Obviously, something is wrong when two divisions with exactly the same operands produce different results, however tiny the difference is. Even worse, as demonstrated by @urbanjost’s experiments, the results become consistent after some seemingly irrelevant changes to the code. At least, this should not be the default behavior. If this behavior is acceptable by default, I cannot conceive how/when to trust the results computed by any code involving division.

The bug has been submitted to Intel, although no attention has been received. Maybe some Intel colleagues on this discourse can file a formal bug report? Many thanks.

19

I would expect a compiler to compute the reciprocal a single time, and then multiply that reciprocal by the array in the array expression. So even if that reciprocal were computed inaccurately, the subsequent multiplications should all be consistent.

I first thought of rounding mode issues. Maybe the array expression is split up into parallel execution segments, and some segments have one rounding mode set and other segments have another rounding mode set.

Another possibility might be inconsistent use of the 80-bit register. Some parallel execution segments might be using the 80-bit register differently than others.

1 Like
20

I have a guess as to what’s happening. you theoretically can get higher throughout out of a CPU by using both multiplication and division so they might have taken the inverse and multiplied by the inverse for some of the columns and divided by the number for others.

1 Like