Ifort (IFORT) 2021.8.0: `1.0E+37 / 1.0E+38 = 0`

[zaikunzhang]

@RonShepard

I think it was useful. There is a range of floating point numbers between 1.0/tiny(1.0) and huge(1.0) where one must be careful with reciprocals and divisions. This is not specific to ifort, this information applies to any compiler. I used ifort and gfortran to demonstrate what was happening because their default settings, but there are almost certainly compiler options for both that would switch the behavior regarding denormals.

I entirely agree with you.

This is why I do not think it is the best strategy to stop the discussion by simply claiming that “to make floating-point calculations work as expected, apply -fp-model = strict”. Logically, such a claim implies that “floating-point calculations are arbitrary / unpredictable without fp-model = strict”. It would be scary if we believe both the validity of the claim and the usefulness of Fortran, because there is almost surely some Fortran code crucial for society running without fp-model = strict.

If you only write code for yourself, you just need to remember fp-model = strict, but how would you make sure that your users remember it, provided that you are developing a piece of software? I do not think there is a way.

Thanks to the insights contained in this discussion, we know precisely how to make the code work as expected without relying on fp-model = strict, and moreover, the insights are applicable to much more general situations. I find this particularly educational — maybe it is me who is ignorant, but I have learned a lot from the discussion, which would not have happened had I taken the comfortable and lazy shortcut of fp-model = strict.

Don’t get me wrong. I am not criticizing fp-model = strict. Indeed, I will remember to apply it whenever dealing with floating-point calculations in production code according to what I have learned here.

Thank you for everything.

[FortranFan]

@zaikunzhang

I hope this discussion is useful or interesting to other users of ifort.

If users pay attention to @greenrongreen 's response in the other thread, this won’t be of much use.

Intel’s reply in that thread suggests with default REAL intrinsic type when optimization comes into play with vectorization, the floating-point operations are generally fraught, particularly with division, whether it is with a normal denominator or a denormal one. And that Intel makes certain recommendations such as use -prec-div.

To further illustrate Intel’s point, consider another case:

   integer, parameter :: WP = selected_real_kind( p=6 )
   real(WP), parameter :: FOO = 0.0999_wp*huge(1.0_wp)
   real(WP), parameter :: BAR = 0.999_wp*huge(1.0_wp)
   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

• Without optimization, default REAL, and default floating-point setting

C:\temp>ifort /standard-semantics /Od 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 =  3.3994206E+37 3.3994206E+37 3.3994206E+37
 After division: x =  0.1000000 0.1000000 0.1000000

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

C:\temp>

• With optimization, default REAL, and default floating-point

C:\temp>ifort /standard-semantics /O2 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 =  3.3994206E+37 3.3994206E+37 3.3994206E+37
 After division: x =  0.1000000 0.1000000 0.1000000

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

C:\temp>

• With optimization, default REAL, and prec-div

C:\temp>ifort /standard-semantics /O2 /Qprec-div 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 =  3.3994206E+37 3.3994206E+37 3.3994206E+37
 After division: x =  0.1000000 0.1000000 0.1000000

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

C:\temp>

• With optimization, double precision, and default floating-point setting

C:\temp>ifort /standard-semantics /O2 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.795895441727453E+307 1.795895441727453E+307 1.795895441727453E+307
 After division: x =  0.100000000000000 0.100000000000000
 0.100000000000000

 Block 2: vector division
 x =  1.795895441727453E+307 1.795895441727453E+307 1.795895441727453E+307
 After division: x =  0.100000000000000 0.100000000000000
 0.100000000000000

C:\temp>

It is the same underlying issue as that other thread, Intel users are not learning anything new here.

[zaikunzhang]

@FortranFan

It is the same underlying issue as that other thread, Intel users are not learning anything new here.

While I am happy to get confirmation that at least someone finds it interesting and useful,
I feel sorry that you find it uninteresting and unuseful but still spend time writing and reading lengthy responses during this season. Thank you for your efforts.

[RonShepard]

This is the relevant point that was exposed in this discussion.

integer, parameter :: WP = selected_real_kind( p=6 )
real(WP), parameter :: BAR = 0.999_wp*huge(1.0_wp)
real(wp), parameter :: big = 1.0_wp / tiny(1.0_wp)
if ( BAR > big ) then
   write(*,*) 'Beware of recipricols and divisons with BAR=', BAR, ' big=', big
endif
end program

When numbers are in this range, >big, then the expressions should be rewritten to avoid the divisions or reciprocals, or the whole problem should be scaled so that such problems do not arise.

[FortranFan]

@RonShepard

When numbers are in this range, >big, then the expressions should be rewritten to avoid the divisions or reciprocals, or the whole problem should be scaled so that such problems do not arise.

Processors have generally gotten much better and for a couple of decades or more, it has generally not been necessary to be so over-analytical with floating-point arithmetic.

Besides, for users of IFORT, it does not explain nor help with the underlying issue, for if it did, a code such as below will throw unexpected behavior regardless of floating-point precision or the optimization setting. But instead it only does so when default REAL is used with Intel’s fast option for floating-point is in effect with optimization being turned on:

   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(5)
   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
      x(4) = x(4) / BAR ; x(5) = x(5) / 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

• Expected behavior per OP when default REAL is used with NO optimization and fp:fast

C:\temp>ifort /standard-semantics /fp:fast /Od q.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:q.exe
-subsystem:console
q.obj

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

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

C:\temp>

• Unexpected behavior per OP when optimization is used with default REAL and fp:fast

C:\temp>ifort /standard-semantics /fp:fast /O2 q.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:q.exe
-subsystem:console
q.obj

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

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

C:\temp>

• Expected behavior per OP if double precision used with optimization and fp:fast option

C:\temp>ifort /standard-semantics /fp:fast /O2 q.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:q.exe
-subsystem:console
q.obj

C:\temp>q.exe
 Block 1: unrolled division
 x =  5.393079404586948E+306 5.393079404586948E+306 5.393079404586948E+306
 5.393079404586948E+306 5.393079404586948E+306
 After division: x =  0.100000000000000 0.100000000000000
 0.100000000000000 0.100000000000000 0.100000000000000

 Block 2: vector division
 x =  5.393079404586948E+306 5.393079404586948E+306 5.393079404586948E+306
 5.393079404586948E+306 5.393079404586948E+306
 After division: x =  0.100000000000000 0.100000000000000
 0.100000000000000 0.100000000000000 0.100000000000000

C:\temp>

And again, with -prec-div, the program behavior per OP will be as expected with all cases including when default REAL is used with fp:fast and optimization.

Bottom-line: The issue, again, has to do with Intel Fortran’s optimization approach with vectorization involving its default REAL type. This violates POLA principle for users such as OP are taken aback to make multiple posts, to write the least.

[urbanjost]

Well, I added tests based on this to subroutines I include for optional execution in some programs that was originally based on the paranoid code by Richard Karpinski, so I have
found this interesting.

And it has
made me reconsider the defaults used with some compilers by fpm(1), and at least until compiler options can be specified in the fpm(1) manifest file it currently only gives two default sets of compiler options. Makes me think there should be a --profile precise option or something similar,
and maybe a suite of test programs that should pass at least the “fpm --profile debug” build for each compiler.

The paranoid procedures derived from Karpinski’s tests are compiled with the same options as the code, and then can optionally be called to test the options you compiled with; which I have found useful in the past; but the ones I have public (in the General Purpose Fortran collection) have not been updated in quite a while and, like the original program, are interactive; not totally automated.

I am wondering whether a public repository capturing some of the points made here with example codes and the great explanations should be developed, including some of the publicly available materials on the perils of floating point ops. Considering how many codes work without issues seems to indicate the defaults of the various compilers are generally fine; and most people do not often work with numbers in the ranges where denormalization and other issues are likely to occur; but when they do they can generate a lot of red herrings.

I am somewhat surprised several of these problem cases went unflagged when I turned on all the traps I could find on a couple of compilers; but there are so many I probably missed some.
From the documentation I thought they would alert on the denormalized values, underflows and overflows, etc… but they appeared to not always catch them with what I tried so far.

If anyone starts a site dedicated to just these issues let me know. This has given me several ideas on things I want to do, but I know I have no time for them at the moment.

Note that the Intel compilers have a feature where you can change the default compiler options globally and via environment variables that we used to use extensively at a previous employer.
Papers were generated based on over two hundred programs and a special version of paranoi was generated that resulted in a preferred set of defaults for the work being conducted, and materials describing why you would change them and the defaults were changed accordingly.
The only one that there was much controversy about was changing the default build to be static.
Last I knew that was a feature unique to ifort, but it may be available with other compilers now as well. Unfortunately, I no longer have access to those materials, but I often wish I still did.

[RonShepard]

Here is a version of your code that incorporates the features of this discussion.

program fp
   implicit none
   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(5), r
   integer :: iexp
   intrinsic :: EXPONENT, SCALE
   print *, 'BIG=', BIG, ' BAR=', BAR
   if ( BAR > BIG ) then
      print *, 'Warning: Beware of divisions and reciprocals with BAR.'
      iexp = 16 - EXPONENT(BAR)
   else
      print *, 'Notice: Divisions and reciprocals are safe with BAR.'
      iexp = 0
   end if
   print *
   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
   x(4) = x(4) / BAR ; x(5) = x(5) / 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
   print *
   blk3: block
   print *, "Block 3: unrolled reciprocal"
   x = FOO
   print *, "x = ", x
   r = 1.0_wp / BAR
   x(1) = x(1) * r ; x(2) = x(2) * r ; x(3) = x(3) * r
   x(4) = x(4) * r ; x(5) = x(5) * r
   print *, "After multiplication x = ", x
   end block blk3
   print *
   blk4: block
   print *, "Block 4: vector reciprocal"
   x = FOO
   print *, "x = ", x
   r = 1.0 / BAR
   x = x * r
   print *, "After multiplication: x = ", x
   end block blk4
   print *
   blk5: block
   print *, 'Block5: safe scaled vector division with iexp=', iexp
   x = FOO
   print *, "x = ", x
   x = SCALE( x, iexp ) / SCALE( BAR, iexp )
   print *, "After division: x = ", x
   end block blk5
   print *
   blk6: block
   print *, 'Block6: safe scaled vector reciprocal with iexp=', iexp
   x = FOO
   print *, "x = ", x
   r = 1.0_wp / SCALE( BAR, iexp )
   x = SCALE( x, iexp ) * r
   print *, "After multiplication: x = ", x
   end block blk6
end program fp

With ifort default options, here is the output.

$ ifort  fp.f90 && a.out
 BIG=  8.5070592E+37  BAR=  1.0208471E+38
 Warning: Beware of divisions and reciprocals with BAR.
 
 Block 1: unrolled division
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block 2: vector division
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.1000000    
 
 Block 3: unrolled reciprocal
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication x =   0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00
 
 Block 4: vector reciprocal
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication: x =   0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.0000000E+00
 
 Block5: safe scaled vector division with iexp=        -111
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block6: safe scaled vector reciprocal with iexp=        -111
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000

As I think is clear, the problems are entirely due to whether the reciprocal+multiplication is used rather than the division when evaluating the expression. When the reciprocal is used, then the result depends on whether the denormals are set to zero. If the operands are scaled as in the last two cases, there are no problems, no surprises, with either expression.

I do not believe this is unique to ifort. I think with the appropriate optimization and compiler options, one can probably get any compiler to produce this same output. I personally also think that this is not a compiler bug. I think the ifort default choices are legitimate, and when these kinds of issues occur due to floating point quirks, it is up to the programmer to fix them. In general, this can be done by rearranging expressions to avoid the overflow/underflow situations, or by scaling the whole problem to avoid such issues.

A developer has the choice of fixing these kinds of problems once and forever in the source code, or of using the right magical combination of compiler options every time the code is compiled. Since it is impractical for a programmer to specify the correct compiler options for every current and future compiler, it is best for the problem to be corrected in the source code. And preferably with sufficient comments and explanations so that other programmers do not mistakenly change the code and reintroduce the problem in the future.

I also think, as a quality of implementation issue, that when such denormal operations occur that the compiler should somehow warn the programmer at run time. This allows the programmer to fix the problem. When these things occur silently, they are a land mine waiting to explode at the worst possible time. Certainly there are many situations where denormals can be set to zero without causing problems, but it should be up to the programmer to make that decision, it should not be done silently by the compiler.

[RonShepard]

Just for comparison here is the gfortran output with default options.

gfortran  fp.f90 && a.out
 BIG=   8.50705917E+37  BAR=   1.02084714E+38
 Warning: Beware of divisions and reciprocals with BAR.

 Block 1: unrolled division
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block 2: vector division
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block 3: unrolled reciprocal
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication x =    9.99999866E-02   9.99999866E-02   9.99999866E-02   9.99999866E-02   9.99999866E-02

 Block 4: vector reciprocal
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication: x =    9.99999866E-02   9.99999866E-02   9.99999866E-02   9.99999866E-02   9.99999866E-02

 Block5: safe scaled vector division with iexp=        -111
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block6: safe scaled vector reciprocal with iexp=        -111
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

The difference, compared to ifort, is that by default gfortran does not set the denormals to zero. Note however, that the unscaled reciprocal expressions result in slightly different answers than the others. The difference is just in the last bit, so it could be a result of the denormal value losing some bits of accuracy, or it could be just normal fp rounding. Notice in the last two outputs, the scaled results happen to have the full accuracy. That is because the mantissa bits are unchanged by the scaling operation and no bits are lost in computing the reciprocal compared to the division. With other values of FOO and BAR, there will be differences in the last bit of the results for the divisions and the scaled reciprocals because of rounding. This is what I would consider expected fp behavior, compared to the large difference observed when the denormals are set to zero (as ifort does).

I just checked the compiler options for gfortran to see if there was a way to set the denormals to zero, and I don’t see anything obvious to try. Also, my version of gfortran does not have the intrinsic ieee_set_underflow() routine, so I can’t see what gfortran does with that either. If someone else knows how to set that option, it would be nice to see the corresponding gfortran output.

[zaikunzhang]

@RonShepard

I just checked the compiler options for gfortran to see if there was a way to set the denormals to zero, and I don’t see anything obvious to try. Also, my version of gfortran does not have the intrinsic ieee_set_underflow() routine, so I can’t see what gfortran does with that either. If someone else knows how to set that option, it would be nice to see the corresponding gfortran output.

It seems that -funsafe-math-optimizations can flush denormals to zero. However, the official documentation of this option is “Allow optimizations that may be give incorrect results for certain IEEE inputs”, without specifying the desired behavior. So I do not know whether this is the best way.

Code:

program test                                                                                            
implicit none                                                                                           
                                                                                                        
write (*, *) 1.0 / 1.0E+38                                                                                  
                                                                                                        
end program test                                                                                                            

Experiement:

$ gfortran --version 
GNU Fortran (Ubuntu 11.3.0-1ubuntu1~22.04) 11.3.0
Copyright (C) 2021 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

$ gfortran test.f90  && ./a.out 

   9.99999935E-39

$ gfortran -funsafe-math-optimizations test.f90  && ./a.out 

   0.00000000  

In contrast, ifort flushes denormals to zero by default, since -ftz is the default:

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

  0.0000000E+00

Nevertheless, ifx does not do the same as ifort.

$ ifx --version && ifx test.f90 && ./a.out 
ifx (IFORT) 2023.0.0 20221201
Copyright (C) 1985-2022 Intel Corporation. All rights reserved.

  1.0000001E-38

[oscardssmith]

Yeah. Intel doesn’t know how to make subnormals fast (unlike everyone else) so they just write compilers that disable them by default.

[RonShepard]

Is it possible to make arithmetic with subnormals (gradual underflow) fast?

[oscardssmith]

On AMD and some Arm chips it’s fast. Intel just missed the memo. See http://www.acsel-lab.com/arithmetic/arith16/papers/ARITH16_Schwarz.pdf for a 2 decade old paper on how to do it.

[RonShepard]

@zaikunzhang

It seems that -funsafe-math-optimizations can flush denormals to zero.

Here is the output for that option:

$ gfortran -funsafe-math-optimizations fp.f90 && a.out
 BIG=   8.50705917E+37  BAR=   1.02084714E+38
 Warning: Beware of divisions and reciprocals with BAR.

 Block 1: unrolled division
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block 2: vector division
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block 3: unrolled reciprocal
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication x =    0.00000000       0.00000000       0.00000000       0.00000000       0.00000000    

 Block 4: vector reciprocal
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication: x =    0.00000000       0.00000000       0.00000000       0.00000000       0.00000000    

 Block5: safe scaled vector division with iexp=        -111
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After division: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

 Block6: safe scaled vector reciprocal with iexp=        -111
 x =    1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37   1.02084708E+37
 After multiplication: x =    9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02   9.99999940E-02

That does make the cases with the unscaled reciprocal flush to zero, the same way that ifort does by default.

[RonShepard]

@oscardssmith

On AMD and some Arm chips it’s fast.

Ah, I misunderstood. When you said “intel” I thought you were talking about ifort and ifx compilers, but you were talking about the hardware.

[RonShepard]

Here is the ifort output with the lines

use, intrinsic :: ieee_arithmetic, only: ieee_set_underflow_mode
...
call ieee_set_underflow_mode(.true.)

added to the above code. This manually turns on the gradual underflow mode for the floating point arithmetic.

$ ifort fp.f90 && a.out
 BIG=  8.5070592E+37  BAR=  1.0208471E+38
 Warning: Beware of divisions and reciprocals with BAR.
 
 Block 1: unrolled division
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block 2: vector division
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.0000000E+00  0.0000000E+00  0.0000000E+00
  0.0000000E+00  0.1000000    
 
 Block 3: unrolled reciprocal
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block 4: vector reciprocal
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block5: safe scaled vector division with iexp=        -111
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After division: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000    
 
 Block6: safe scaled vector reciprocal with iexp=        -111
 x =   1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37  1.0208471E+37
 After multiplication: x =   0.1000000      0.1000000      0.1000000    
  0.1000000      0.1000000

This is still with the default ifort options. Note that the vector division case is still inconsistent, while the cases that use the denormal reciprocal are now correct. Presumably that is because the AVX hardware that computes the first four terms in the vector expression is unaffected by the gradual underflow setting.

It is also interesting that if I switch off the gradual underflow and compile with gfortran, the output remains unchanged. Both these ifort and gfortran calculations are run with MacOS on intel hardware.

[edit]
I have looked at this behavior a bit more. ieee_get_underflow_mode(gradual) on gfortran shows that the mode is being set correctly, but the output from the reciprocal cases does not reflect that setting. The same output results from either .true. or .false.. These calculations were done with MacOS and gfortran 11.3 (I think the same thing occurs with 12.x).

For ifort, the underflow mode is also set correctly, and the output reflects that setting as expected. Namely, for the reciprocal expressions with gradual =.false., the output values are zero, indicating that the denormals are set to zero prior to the multiplication step.

For ieee arithmetic, I think huge(1.0_wp)*tiny(1.0_wp) is 4.0_wp for all real kinds. Or maybe it is nearest(4.0_wp,-1.0_wp), I’m uncertain what that last bit is really supposed to be. In any case, that means that the large floating point values with the largest four exponents result in denormal reciprocal values, but only the least significant 1 or 2 bits might be incorrect, so there should be only a minor loss of precision in these division expressions when gradual==.true..