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.