How to avoid underflow when multiplying real numbers?

Question

I would like to know if there is a proper way to multiply or square float (or double) numbers together without having underflow error when I compile my Fortran code like

gfortran -ffpe-trap=invalid,zero,overflow,underflow ...

I know that the underflow option is not always a good idea, but I wonder if it is possible to do multiplication with this option. In fact, in the following example I know that underflow may occur but maybe I'm not aware of other case in my codes. This is why I would like to keep this option, if possible.

Here is a example where I compute a vector u for every x,y indices of a matrix; the 2 values composing theses vectors are between 0 and 1. Then I compute the square of its norm.

So very logical, I will have values underflowing because of this square operation. Since, these very small values can be considered as zero for me. Is there a way to not underflow better than using an if comparison?

implicit none
double  :: u(100,100,2), uSqr(100,100)
integer :: x,y

DO x= 1, 100
    DO y = 1, 100
                CALL Poisin( u(x,y,:), x, y )
    ENDDO
ENDDO

uSqr = u(:,:,1)*u(:,:,1) + u(:,:,2) * u(:,:,2) ! where comes the underflow errors

Answer 1

You have an answer which looks at a specific way of avoiding undue underflow in certain circumstances. This uses the hypot function. This is partly an answer: if you want to avoid underflow there may be a way to rewrite an algorithm to avoid it.

For more general cases (such as in this question) where fine control of exception flags is wanted that won't be suitable. However, compilers often offer interfaces to exception-handling routines.

One portable way of doing this is using the IEEE facility of Fortran 2003. [If using gfortran you'll need at least version 5.0, but there are similar compiler-specific ways available.]

Fortran defines IEEE exceptions and flags. A flag can be quiet or signaling. What you want is for the part where underflow is not a useful diagnostic to not affect the underflow flag status after that computation.

The flag is known as IEEE_UNDERFLOW. We can query and set its status with the subroutines calls IEEE_GET_FLAG(IEEE_UNDERFLOW, value) and IEEE_SET_FLAG(IEEE_UNDERFLOW, value). If we're expecting, but don't care about, an underflow, we'll also want to ensure the exception is non-halting. The subroutine IEEE_SET_HALTING_MODE(IEEE_UNDERFLOW, value) controls this mode.

So, an annotated example.

  use, intrinsic :: ieee_arithmetic, only : IEEE_SELECTED_REAL_KIND
  use, intrinsic :: ieee_exceptions

  implicit none

  ! We want an IEEE kind, but this doesn't ensure support for underflow control
  integer, parameter :: rk=IEEE_SELECTED_REAL_KIND(6, 70)

  ! State preservation/restoration
  logical should_halt, was_flagged

  real(rk) x

  ! Get the original halting mode and signal state 
  call ieee_get_halting_mode(IEEE_UNDERFLOW, should_halt)
  call ieee_get_flag(IEEE_UNDERFLOW, was_flagged)

  ! Ensure we aren't going to halt on underflow
  call ieee_set_halting_mode(IEEE_UNDERFLOW, .FALSE.)

  ! The irrelevant computation   
  x=TINY(x)
  x=x**2

  ! And restore our old state
  call ieee_set_halting_mode(IEEE_UNDERFLOW, should_halt)
  call ieee_set_flag(IEEE_UNDERFLOW, was_flagged)

end program