Checking for Additive Underflow

The code shown below is the emulation function for the addition of two emulated 4-byte real numbers. As in the emulated multiplication, the precision of the emulated numbers is adjusted by the routine experiment_4. After the result has been computed, the check for underflow is made by the routine check_underflow_r4_r4.

ELEMENTAL FUNCTION add_em_real_k4_em_real_k4(a1mr4,a2mr4) IMPLICIT NONE TYPE (em_real_k4) :: add_em_real_k4_em_real_k4 TYPE (em_real_k4),INTENT(IN) :: a1mr4 TYPE (em_real_k4),INTENT(IN) :: a2mr4 TYPE (em_real_k4) :: t1mr4 TYPE (em_real_k4) :: t2mr4 t1mr4 = experiment_4(a1mr4) t2mr4 = experiment_4(a2mr4) add_em_real_k4_em_real_k4%value = t1mr4%value + t2mr4%value add_em_real_k4_em_real_k4 = experiment_4(add_em_real_k4_em_real_k4) CALL check_underflow_r4_r4(add_em_real_k4_em_real_k4%value, & a1mr4%value,a2mr4%value) END FUNCTION add_em_real_k4_em_real_k4

The routine check_underflow_r4_r4 checks for additive underflow in the addition or subtraction of two 4-byte numbers. This kind of underflow occurs when the magnitude of the two input numbers differs by so much that the mantissa of the number with the smaller magnitude is severely truncated in the addition or subtraction.

In this routine, a and b are the two numbers added or subtracted, and r is the result. Two criteria are applied in the check. If the magnitude of a is greater than that of b, then no check is made if the magnitude of b is less than a criterion value of check_underflow_delta. This prevents false positives from being reported when the smaller magnitude is very close to zero. We then compute the difference between the result of the addition or subtraction and the value of b. If this differs by a criterion proportion from the magnitude of b then a large part of b has been lost in the computation and a report is made.

Note that the report is made by setting values in shared memory which are detected by the trace routines described earlier. The underflow check cannot write a message directly because elemental sub-programs may not contain I/O statements.

ELEMENTAL SUBROUTINE check_underflow_r4_r4(r,a,b) IMPLICIT NONE REAL(KIND=kr4),INTENT(IN) :: r REAL(KIND=kr4),INTENT(IN) :: a REAL(KIND=kr4),INTENT(IN) :: b REAL (KIND=kr8) :: ta,tb,tmax,tmin,tr,td,te ta = ABS(a) tb = ABS(b) IF (ta > tb) THEN tmax = ta tmin = tb ELSE tmax = tb tmin = ta ENDIF tr = ABS(r) td = ABS(tr-tmax) IF (tmin > check_underflow_delta) THEN te = ABS(1.0D0-td/tmin) IF (te > check_underflow_crit) THEN CALL signal_underflow(te,tr,tmax,tmin) ENDIF ENDIF END SUBROUTINE check_underflow_r4_r4