Exception handling in Fortran
John Reid
14 October 1995
I tried hard to get an exception handling feature into Fortran during
1993 and 1994. Although it was accepted by the ISO committee (WG5) in
August 1994, there were detailed criticisms by members of the ANSI
committee (X3J3) which led X3J3 to recommend the abandonment of the
feature for Fortran 95. There was a worry that addressing the
criticisms might lead to a delay to the whole standard.
This recommendation was endorsed by WG5, but at its April 1995 meeting in
Tokyo it decided that handling floating-point exceptions was too
important to leave until Fortran 2000. It therefore decided to
establish a development body to create a "Type 2 Technical Report". The
intention is to finalize this soon (to be ready for formal balloting by
April 1996). It will permit vendors to implement the feature as an
extension of Fortran 95, confident that it will be part of Fortran
2000, unless experience in their implementation and use demand changes.
It is a kind of beta-test facility for a new language feature.
I have agreed to act as project editor and have been asked to provide a
draft based on last year's work and an alternative proposal using
intrinsic procedures. The intrinsics approach has always been advocated
by Kahan, but I have not pursued it in the past because everyone has
found the enable approach much more elegant. However, it is difficult
to understand the fine details of the enable approach and different
people have very different expectations of it. The intrinsics are easy
to understand and my current thinking is that this is far more likely
to be acceptable to people.
You have seen the enable proposal before, so I will confine myself to
the intrinsics proposal, whose design is due to Keith Biermann of SUN.
It addresses more than just the exceptions, and I hope that WG 2.5 will
support the approach.
The idea is to have an intrinsic module containing the types
1. IEEE_FLAG, with the constants INVALID, UNDERFLOW, OVERFLOW,
INEXACT, DIVIDE_BY_ZERO, COMMON. (COMMON combines INVALID,
OVERFLOW, and DIVIDE_BY_ZERO.)
2. IEEE_RND_VALUE, with the constants NEAREST, TO_ZERO, TO_INFINITY.
3. IEEE_FSR_VALUE, for saving the current state of the
floating point environment (usually the floating point
status register).
There will be the following procedures:
A. Inquiry functions:
IEEE_SUPPORT_ALL() Inquire if processor supports all the IEEE facilities
defined in this standard.
IEEE_SUPPORT_NAN(X) Inquire if processor supports the IEEE Not-A-Number
facility for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_INF(X) Inquire if processor supports the IEEE infinity facility
for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_FLAG(flag) Inquire if the processor supports flag.
IEEE_SUPPORT_DENORMAL(X) Inquire if the processor supports IEEE gradual
underflow for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_SQRT(X) Inquire if the processor supports IEEE square root
for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_HALTING() Inquire if the processor supports the ability
to control during program execution whether to abort or continue
execution after IEEE exceptions.
IEEE_SUPPORT_ROUNDING (RND_VALUE) Inquire if processor supports aZERO
particular IEEE rounding mode. RND_VALUE is of type IEEE_RND_VALUE,
with value NEAREST, TO_ZERO, or TO_INFINITY.
IEEE_SUPPORT_FSR() Inquire whether the processor supports getting and
setting of the set of flags that define the current state of the
floating point environment (usually in the floating point status
register).
IEEE_DATATYPE (X) Inquire if the processor supports IEEE arithmetic
for reals of the same kind type parameter as the argument.
B. Elemental functions:
IEEE_SQRT(X) Square root as defined in the IEEE standard.
IEEE_IS_NAN(X) Determine if value is IEEE Not-a-Number.
IEEE_IS_INF(X) Determine if value is IEEE Infinity.
IEEE_IS_VALID(X) Determine if value is neither an Infinity nor a NaN.
IEEE_IS_DENORMAL(X) Determine if value is IEEE denormalized.
IEEE_COPYSIGN(X,Y) IEEE copysign function.
IEEE_NEXTAFTER(X,Y) Returns the next representable neighbor of X in the
direction toward Y.
IEEE_LOGB(X) Unbiased exponent in the IEEE floating point format.
IEEE_SCALB (X,I) Returns X * 2**I.
IEEE_INFINITY(X) Generate IEEE infinity with the same sign as X.
IEEE_NAN(X) Generate IEEE Not-a-Number.
IEEE_UNORDERED(X,Y) IEEE unordered function. True if X or Y is NaN
and false otherwise.
C. Scalar functions:
IEEE_FLAG_GET(FLAG) Get an IEEE flag (also known as exception or
sticky bit), where FLAG is of type IEEE_FLAG and value INVALID,
UNDERFLOW, OVERFLOW, INEXACT, or DIVIDE_BY_ZERO.
IEEE_GET_FSR() Save the current values of the set of flags that define
the current state of the floating point environment (usually the
floating point status register). The result is of type
IEEE_FSR_VALUE.
IEEE_GET_ROUNDING_MODE() Save the current IEEE rounding mode. The
result is of type IEEE_RND_VALUE.
D. Subroutines:
IEEE_FLAG_SET(flag-list) Set flags listed.
IEEE_FLAGS_CLEAR(flag-list) Clear flags listed.
IEEE_SET_FSR(FSR_VALUE) Restore the values of the set of flags that
define the current state of the floating point environment (usually
the floating point status register). FSR_VALUE is of type
IEEE_FSR_VALUE and has been set by a call of IEEE_GET_FSR.
IEEE_SET_ROUNDING_MODE(RND_VALUE) Set the current IEEE rounding mode.
RND_VALUE is of type IEEE_RND_VALUE, with value NEAREST, TO_ZERO, or
TO_INFINITY.
IEEE_HALT(FLAG, HALTING) Controls continuation or halting on exceptions.
FLAG is of type IEEE_FLAG and value INVALID, UNDERFLOW, OVERFLOW,
INEXACT, or DIVIDE_BY_ZERO.
I think the minimum requirement is for the support of OVERFLOW and
DIVIDE_BY_ZERO. All the inquiry functions need to be there (but could all
return false) and the procedures for getting, setting, and clearing the
flags need to be there.
Something needs to be said about what happens in PURE procedures. My
current thought is that if a processor spawns tasks in a FORALL:
a. the processor should store the flags before executing the FORALL and
copy them to each processor;
b. at each join (when a processor is done) the stored flags should
be OR'd with the processor values.
Example 1:
USE IEEE_ARITHMETIC
TYPE(IEEE_FRS_VALUE) FSR_VALUE
FSR_VALUE = IEEE_GET_FSR()
CALL IEEE_FLAGS_CLEAR (COMMON)
! First try the "fast" algorithm for inverting a matrix:
MATRIX1 = FAST_INV (MATRIX)
! MATRIX is not altered during execution of FAST_INV.
IF (IEEE_FLAG_GET(COMMON)) THEN
! "Fast" algorithm failed; try "slow" one:
CALL IEEE_FLAGS_CLEAR (OVERFLOW, INVALID, DIVIDE_BY_ZERO)
MATRIX1 = SLOW_INV (MATRIX)
IF (IEEE_FLAG_GET(COMMON) ) THEN
WRITE (*, *) 'Cannot invert matrix'
STOP
END IF
CALL IEEE_SET_FSR(FSR_VALUE)
In this example, the function FAST_INV may cause a condition to
signal. If it does, another try is made with SLOW_INV. If this still
fails, a message is printed and the program stops. Note, also, that it
is important to clear the flags before the second try.
Example 2 (from Tom Hull, modified to use these features):
REAL FUNCTION HYPOT(X, Y)
! In rare circumstances this may lead to the setting of the OVERFLOW flag
USE IEEE_ARITHMETIC
REAL X, Y
REAL SCALED_X, SCALED_Y, SCALED_RESULT
LOGICAL OLD_OVERFLOW, OLD_UNDERFLOW
INTRINSIC SQRT, ABS, EXPONENT, MAX, DIGITS, SCALE
! Store the old flags and clear them
OLD_OVERFLOW = IEEE_FLAG_GET(OVERFLOW)
OLD_UNDERFLOW = IEEE_FLAG_GET(UNDERFLOW)
CALL IEEE_FLAGS_CLEAR(UNDERFLOW, OVERFLOW)
! Try a fast algorithm first
HYPOT = SQRT( X**2 + Y**2 )
IF (IEEE_FLAG_GET(UNDERFLOW) .OR. IEEE_FLAG_GET(OVERFLOW) ) THEN
CALL IEEE_FLAGS_CLEAR(UNDERFLOW, OVERFLOW)
IF ( X==0.0 .OR. Y==0.0 ) THEN
HYPOT = ABS(X) + ABS(Y)
ELSE IF ( 2*ABS(EXPONENT(X)-EXPONENT(Y)) > DIGITS(X)+1 ) THEN
HYPOT = MAX( ABS(X), ABS(Y) )! one of X and Y can be ignored
ELSE ! scale so that ABS(X) is near 1
SCALED_X = SCALE( X, -EXPONENT(X) )
SCALED_Y = SCALE( Y, -EXPONENT(X) )
SCALED_RESULT = SQRT( SCALED_X**2 + SCALED_Y**2 )
HYPOT = SCALE( SCALED_RESULT, EXPONENT(X) ) ! may cause overflow
END IF
END IF
IF(OLD_OVERFLOW) CALL IEEE_FLAG_SET(OVERFLOW)
IF(OLD_UNDERFLOW) CALL IEEE_FLAG_SET(UNDERFLOW)
END FUNCTION HYPOT
An attempt is made to evaluate this function directly in the fastest
possible way. (Note that with hardware support, exception checking is
very efficient; without language facilities, reliable code would
require programming checks that slow the computation significantly.)
The fast algorithm will work almost every time, but if an exception
occurs during this fast computation, a safe but slower way evaluates
the function. This slower evaluation may involve scaling and
unscaling, and in (very rare) extreme cases this unscaling can cause
overflow (after all, the true result might overflow if X and Y are both
near the overflow limit). If the overflow or underflow flag is set on
entry, it is reset on return, so that earlier exceptions are not
lost.