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

To be clear, do you want to avoid underflow, or do you want to avoid having underflow being an error condition? — francescalus, Apr 11 '17 at 17:46
Can't you set the `-ffpe-trap` option on a per-file basis? If so, maybe you could separate out the part cited above into one source file and omit `underflow` from the`-ffpe-trap` list to compile this file while keeping the original exception trap list for the other source files. I assume that you want to allow the program to go on despite the underflow of `uSqr` values by regarding them as zeros. — norio, Apr 11 '17 at 18:08
I must second @francescalus's comment. Do you want to keep the computation and just not exit with an error? Why do you want or need to trap `underflow`? Seriously consider not trapping `underflow`. In my opinion it is pretty useless or even harmful. — Vladimir F Героям слава, Apr 11 '17 at 21:08
@norio This is what I do know. But with this solution I don't know if their is some other `underflow` error in this file. — R. N, Apr 12 '17 at 17:39
@VladimirF and @francescalus I want to trap `underflow` to be sure their is no `underflow` error in the rest of my code. So I want to know if their is a better to write my example that avoids `underflow` problems. Maybe my question is stupid and the answer is because I do square operations `underflow` and `overflow` might occurs so no their is no other way except if I use a conditional test to threshold the 'small' and 'big' values. — R. N, Apr 12 '17 at 17:39

score 3 · Accepted Answer · answered Apr 12 '17 at 17:43

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

your answer seem very appropriate. I will try at home this weekend and I will validate your answer. — R. N, Apr 14 '17 at 10:09

score 1 · Answer 2 · answered Apr 12 '17 at 08:32

1

If it's just to avoid underflow maybe you want to use hypot or do you really need the square of hypotenuse?

The implementation of hypot should avoid overflow/underflow problems of sqrt(x**2+y**2).

answered Apr 12 '17 at 08:32

aka.nice

9,100
1
28
40

You are right, I always considered it only for overflow problems, but it should avoid underflow in the same way. – Vladimir F Героям слава Apr 12 '17 at 10:38
Thanks for the `hypot` tips, but I really need the square of the hypotenuse. – R. N Apr 12 '17 at 17:29
@Romain So compute the square of hypot, where is the problem... If the underflow happens during that, it is unavoidable. – Vladimir F Героям слава Apr 12 '17 at 17:41

How to avoid underflow when multiplying real numbers?

2 Answers2