Flush-to-zero in gfortran

2019-06-27 14:50发布

问题:

Is there a way to force flush-to-zero of underflows in gfortran?

I can't believe this is the first time someone has asked this, but I couldn't find anything on it anywhere. Mea culpa if this is a duplicate.

回答1:

You can accomplish this with recent versions of gfortran that support the Fortran 2003 IEEE modules. The standard defines two underflow modes -- gradual and abrupt. Abrupt is the one you want which sets underflow to 0 and signals the underflow floating point exception. You can test for support of controlling the underflow mode with the function ieee_support_underflow_control(X) which tests for underflow control for the kind of real X is and returns a logical true if it is supported. If supported, you can then call ieee_set_underflow_mode(.false.) to set abrupt underflow mode.

Below is a test program you can use to test underflow control support for the default real kind:

program test
  use, intrinsic :: ieee_arithmetic
  use, intrinsic :: iso_fortran_env, only: compiler_version, compiler_options
  implicit none
  logical :: underflow_support, gradual, underflow
  real :: fptest
  integer :: i

  print '(4a)',  'This file was compiled by ', &
       compiler_version(), ' using the options ', &
       compiler_options()
  fptest = 0.0
  underflow_support = ieee_support_underflow_control(fptest)
  if (underflow_support) then
     print *,'Underflow control supported for the default real kind'
  else
     stop 'no underflow control support'
  end if

  call ieee_set_underflow_mode(.false.)
  call ieee_get_underflow_mode(gradual)
  if (.not.gradual) then 
     print *,'Able to set abrupt underflow mode'
  else
     stop 'error setting underflow mode'
  end if

  fptest = 2e-36
  do i=1,50 ! 50 iterations max
     fptest = fptest * 0.5
     print '(e15.10)',fptest
     call ieee_get_flag(ieee_underflow,underflow)
     if (underflow) print *,'Underflow exception signaling'
     if (fptest == 0.0) exit
  end do

end program test

Using gfortran version 5.2.0, this program outputs:

This file was compiled by GCC version 5.2.0 using the options -mtune=generic -march=x86-64 -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans
 Underflow control supported for the default real kind
 Able to set abrubpt underflow mode
.1000000036E-35
.5000000180E-36
.2500000090E-36
.1250000045E-36
.6250000225E-37
.3125000112E-37
.1562500056E-37
.0000000000E+00
 Underflow exception signaling

The compiler option flags -fno-unsafe-math-optimizations -frounding-math -fsignaling-nans are suggested by the gfortran 5.2 documentation to be used anytime the IEEE modules are used to ensure adherence to the standard.



回答2:

A lazy way to "flush to zero" is to use to use gfortran's -funsafe-math-optimizations to:

Allow math optimizations that may violate IEEE or ISO standards

or in other words:

This mode enables optimizations that allow arbitrary reassociations and transformations with no accuracy guarantees. It also does not try to preserve the sign of zeros.

For example, small.f:

      program test
        real r
        r=1e-40
        print *,'r on next line'
        print *,r
      end program

Without any flags, a non-zero denormal (small) number is shown, without error:

$ gfortran -g small.f
$ ./a.out
 r on next line
   9.99994610E-41

Trap the denormalized number, which crashes when attempting to print the value:

$ gfortran -g -ffpe-trap=denorm small.f
$ ./a.out
 r on next line

Program received signal SIGFPE: Floating-point exception - erroneous arithmetic operation.

Backtrace for this error:
#0  0x2aaaab05c26f in ???
#1  0x2aaaaac61aed in get_float_string
        at ../../../libgfortran/io/write_float.def:1064
#2  0x2aaaaac6423d in list_formatted_write_scalar
        at ../../../libgfortran/io/write.c:1889
#3  0x4008f1 in test
        at /path/to/small.f:5
#4  0x400941 in main
        at /path/to/small.f:6
Floating point exception

And added flag that flushes it to zero:

$ gfortran -g -ffpe-trap=denorm -funsafe-math-optimizations small.f
$ ./a.out
 r on next line
   0.00000000