Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Flush-to-zero in gfortran

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.

like image 686
bob.sacamento Avatar asked Sep 29 '15 18:09

bob.sacamento


2 Answers

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.

like image 69
casey Avatar answered Nov 12 '22 01:11

casey


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
like image 29
Mike T Avatar answered Nov 12 '22 00:11

Mike T