From: John Reid
Subject: Procedures technical report
Date: 1 November 1995
0. NOTES
This is a draft Technical Report using procedures in a module for
exception handling. It is based on Keith Bierman's message of 7
September. For IEEE hardware, it supports exceptions and other aspects
of IEEE arithmetic. Non-IEEE processors are required to support only
OVERFLOW and DIVIDE_BY_ZERO, and may support more. We could reduce
this proposal to one that supports exceptions only. This would involve
removal of the data type IEEE_RND_VALUE and the reduction of the set of
procedures to IEEE_GET_FLAG(FLAG), IEEE_CLEAR_FLAGS, IEEE_SET_FLAGS,
IEEE_GET_FSR, IEEE_HALT, and IEEE_SET_FSR.
More work is needed yet on the interaction of these features with
PURE and FORALL.
1. RATIONALE
Exception handling is required for the development of robust and
efficient numerical software. In particular, it is necessary in order
to be able to write portable scientific libraries. In numerical
Fortran programming, current practice is to employ whatever exception
handling mechanisms are provided by the system/vendor. This clearly
inhibits the production of fully portable numerical libraries and
programs. It is particularly frustrating now that IEEE arithmetic is so
widely used, since built into it are the five conditions: overflow,
invalid, divide_by_zero, underflow, and inexact. Our aim to to provide
support for these conditions.
This proposal involves a standard module, IEEE_ARITHMETIC, containing
three derived types, some named constants of these types, and a
collection of simple procedures. They allow the flags to be tested,
cleared, set, saved, or restored. To facilitate maximum performance,
each of the proposed functions does very little processing of
arguments. In many cases, a processor may generate only a few inline
machine code instructions rather than library calls.
In order to allow for the maximum number of processors to provide the
maximum value to users, we do NOT require complete IEEE
conformance. What a vendor must do is to provide the module, and
support the overflow and divide-by-zero flags and the basic inquiry
functions. A user must utilize the inquiry functions to determine if
they can count on specific features of the IEEE standard, before using
other inquiry or value setting procedures.
Note that a processor ought not implement these as "macros", as IEEE
conformance is often controlled by compiler switches. A processor
which offers a switch to turn off a facility should adjust the values
returned for these inquiries. For example, a processor which allows
gradual underflow to be turned off (replaced with flush to zero)
should return .FALSE. for IEEE_SUPPORT_DENORMAL(X) when a source file
is processed with that option on. Naturally it should return
.TRUE. when that option is not in effect.
The most important use of a floating-point exception handling facility
is to make possible the development of much more efficient software
than is otherwise possible. The following "hypotenuse" function,
sqrt(x**2 + y**2), illustrates the use of the facility in developing
efficient software.
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_GET_FLAG(OVERFLOW)
OLD_UNDERFLOW = IEEE_GET_FLAG(UNDERFLOW)
CALL IEEE_CLEAR_FLAGS(UNDERFLOW, OVERFLOW)
! Try a fast algorithm first
HYPOT = SQRT( X**2 + Y**2 )
IF (IEEE_GET_FLAG(UNDERFLOW) .OR. IEEE_GET_FLAG(OVERFLOW) ) THEN
CALL IEEE_CLEAR_FLAGS(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_SET_FLAGS(OVERFLOW)
IF(OLD_UNDERFLOW) CALL IEEE_SET_FLAGS(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.
Can all this be accomplished without the help of an exception handling
facility? Yes, it can - in fact, the alternative code can do the job,
but of course it is much less efficient. That's the point. The HYPOT
function is special, in this respect, in that the normal and
alternative codes try to accomplish the same task. This is not always
the case. In fact, it very often happens that the alternative code
concentrates on handling the exceptional cases and is not able to
handle all of the non-exceptional cases. When this happens, a program
which cannot take advantage of hardware flags could have a structure
like the following:
if ( in the first exceptional region ) then
handle this case
else if ( in the second exceptional region ) then
handle this case
:
else
execute the normal code
end
But this is not only inefficient, it also "inverts" the logic of the
computation. For other examples, see Hull, Fairgrieve and Tang (1994)
and Demmel and Li (1994).
The HYPOT algorithm is used in example 4 of section 15.10. The
code for the HYPOT function can be generalized in an obvious
way to compute the Euclidean Norm, sqrt( x(1)**2 + x(2)**2 + ... +
x(n)**2 ), of an n-vector; the generalization of the alternative code is
not so obvious (though straightforward) and will be much slower
relative to the normal code than is the case with the HYPOT function.
Of the five IEEE conditions, underflow and inexact are obviously milder
than the other three (inexact is particularly mild and signals almost
all the time in typical codes). We therefore group the other three
together as USUAL.
There is a need for further intrinsic conditions in connection with
reliable computation. Examples are
a. INSUFFICIENT_STORAGE for when the processor is unable to find
sufficient storage to continue execution.
b. INTEGER_OVERFLOW and INTEGER_DIVIDE_BY_ZERO for when an
intrinsic integer operation has a very large result or
has a zero denominator.
c. INTRINSIC for when an intrinsic procedure has been unsuccessful.
d. SYSTEM_ERROR for when a system error occurs.
This proposal has been designed to allow such enhancements in the future.
References
Demmel, J.W. and Li, X. (1994). Faster Numerical Algorithms via
Exception Handling. IEEE Transactions on Computers, 43,
no. 8, 983-992.
Hull, T.E., Fairgrieve, T.F., and Tang, T.P.T. (1994).
Implementing complex elementary functions using exception handling.
ACM Trans. Math. Software 20, 215-244.
2. TECHNICAL SPECIFICATION
This proposal involves a standard module, IEEE_ARITHMETIC, containing
three derived types:
1. IEEE_FLAG, with the constants INVALID, OVERFLOW, DIVIDE_BY_ZERO,
USUAL, UNDERFLOW, INEXACT, and ALL.
2. IEEE_RND_VALUE, with the constants NEAREST, TO_ZERO, UP, DOWN, and
OTHER
3. IEEE_FSR_VALUE, for saving the current floating point status.
and a collection of simple procedures.
There must be flags that are initially clear and that are set when an
exception occurs. The value of a flag is determined by the function
IEEE_GET_FLAG(flag)
Flags may be cleared by the subroutine
IEEE_CLEAR_FLAGS(flag-list)
and may be set by the subroutine
IEEE_SET_FLAGS(flag-list)
The minimum requirement is for the support of OVERFLOW and DIVIDE_BY_ZERO.
Whether the other exceptions are supported may be determined by the
inquiry function
IEEE_SUPPORT_FLAG(FLAG)
The extent of further support for the IEEE standard may be determined
by the inquiry functions:
IEEE_DATATYPE (X) Inquire if the processor supports IEEE arithmetic
for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_ALL() Inquire if processor supports all the IEEE facilities
defined in this standard for all kinds of reals.
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_FLAG(FLAG) Inquire if the processor supports an exception
or set of exceptions.
IEEE_SUPPORT_HALTING() Inquire if the processor supports the ability
to control during program execution whether to abort or continue
execution after an exception.
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_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_ROUNDING (RND_VALUE) Inquire if processor supports
changing to a particular rounding mode for IEEE kinds of reals.
IEEE_SUPPORT_SQRT(X) Inquire if the processor supports IEEE square root
for reals of the same kind type parameter as the argument.
Implementations with the appropriate extent of IEEE support provide:
A. Elemental functions:
IEEE_COPYSIGN(X,Y) IEEE copysign function.
IEEE_INFINITY(X) Generate IEEE infinity with the same sign as X.
IEEE_IS_DENORMAL(X) Determine if value is IEEE denormalized.
IEEE_IS_INF(X) Determine if value is IEEE Infinity.
IEEE_IS_NAN(X) Determine if value is IEEE Not-a-Number.
IEEE_IS_NORMAL(X) Whether a value is normal, that is, neither an
Infinity, a NaN, nor denormalized.
IEEE_LOGB(X) Unbiased exponent in the IEEE floating point format.
IEEE_NAN(X) Generate IEEE Not-a-Number.
IEEE_NEXTAFTER(X,Y) Returns the next representable neighbor of X in the
direction toward Y.
IEEE_SCALB (X,I) Returns X * 2**I.
IEEE_SQRT(X) Square root as defined in the IEEE standard.
IEEE_UNORDERED(X,Y) IEEE unordered function. True if X or Y is a NaN
and false otherwise.
B. Scalar functions:
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.
C. Subroutines:
IEEE_CLEAR_FLAGS (FLAG1 [, FLAG2, ...]) Clear one or more exception flags.
IEEE_HALT(FLAG,HALTING) Controls continuation or halting on exceptions.
FLAG is of type IEEE_FLAG and value INVALID, OVERFLOW, DIVIDE_BY_ZERO,
USUAL, UNDERFLOW, INEXACT, or ALL.
IEEE_SET_FLAGS(FLAG1 [, FLAG2, ...]) Set one or more exception flags.
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,
UP, or DOWN.
Note: An IEEE processor is required to support single precision, and
an extended precision. This is usually REAL and DOUBLE PRECISION but
it might not always be. In addition, some vendors support a "REAL*16"
which does not follow the "logical extension" to the IEEE 754 standard
(e.g. IBM's RS6K "doubled double" implementation of REAL*16). Users
should be able to determine if a particular datatype is going to
behave as "one would expect". But implementors should not be
constrained to supply only datatypes which conform to IEEE 754 (for
performance, or compatibility with older processors or eventual
improvements in state of the art in floating point implementation).
If a processor supports a particular IEEE exception, the associated
flag is always clear.
3. EDITS TO THE DRAFT STANDARD (N1122)
130/11. Add: 'If any exception (15) is set, the processor shall issue a
warning on the unit identified by * in a WRITE statement, indicating
which exceptions are set.'.
288+. Add
<<15. Intrinsic module for support of exceptions and IEEE arithmetic>>
The intrinsic module IEEE_ARITHMETIC provides support for the floating
point exceptions associated with overflow and divide by zero for real
or complex data with any kind parameter. If the type and kind parameter
implements IEEE arithmetic (IEEE 754-1985), the module provides more
extensive support, including the exceptions invalid, underflow, and
inexact. There are inquiry functions that permit the extent of support
of IEEE 754-1985 to be determined.
<<15.1 Derived data types defined in the module>>
The module contains three data types, whose components are private.
No operation is defined for any of them and only intrinsic assignment
is available for them. They are:
(i) TYPE (IEEE_FLAG) is used to identify exception flags or sets
of flags. The constants INVALID, OVERFLOW, DIVIDE_BY_ZERO,
USUAL, UNDERFLOW, INEXACT, and ALL are defined for it.
(ii) TYPE (IEEE_RND_VALUE) is used to identify rounding modes.
The constants NEAREST, TO_ZERO, UP, DOWN, and OTHER are defined
for it.
(iii) TYPE (IEEE_FSR_VALUE) is used for saving the current
floating point status. No constants are defined for it.
<<15.2 The exceptions>>
The exceptions are:
OVERFLOW
This exception occurs when the result for an intrinsic real operation
has a very large processor-dependent absolute value, or the real or
imaginary part of the result for an intrinsic complex operation has a
very large processor-dependent absolute value.
DIVIDE_BY_ZERO
This exception occurs when a real or complex division has a nonzero
numerator and a zero denominator.
INVALID
This exception occurs when a real or complex operation is invalid.
UNDERFLOW
This exception occurs when the result for an intrinsic real operation
has a very small processor-dependent absolute value, or the real or
imaginary part of the result for an intrinsic complex operation has a
very small processor-dependent absolute value.
INEXACT
This exception occurs when the result of a real or complex operation
is not exact.
Each exception has a flag whose value is either clear or set. The value
may be determined by the function IEEE_GET_FLAG. Its initial value is
clear and it becomes set when the associated exception occurs. It may
also be set by the subroutine IEEE_SET_FLAGS or the subroutine
IEEE_SET_FSR. Once set, it remains set unless cleared by an invocation
of the subroutine IEEE_CLEAR_FLAGS or the subroutine IEEE_SET_FSR. If
any exception is set when the program terminates, the processor shall
issue a warning on the unit identified by * in a WRITE statement,
indicating which conditions are set.
If a flag is set on entry to a pure procedure, the flag shall not
be cleared unless it is reset before return.
In a sequence of statements that contains no exception handling
statements, if the execution of a process would cause an exception to
be set but after execution of the sequence no value of a variable
depends on the process, whether the exception is set is processor
dependent. For example, when Y has the value zero, whether the code
X = 1.0/Y
X = 3.0
sets DIVIDE_BY_ZERO is processor dependent.
An exception must not set if this could arise only during
execution of a process further to those required or permitted by the
standard. For example, the statement
IF (F(X)>0.) Y = 1.0/Z
must not set DIVIDE_BY_ZERO when both F(X) and Z are zero and the
statement
WHERE(A>0.) A = 1.0/A
must not set DIVIDE_BY_ZERO. On the other hand, when
X has the value 1.0 and Y has the value 0.0, the expression
X>0.00001 .OR. X/Y>0.00001
is permitted to cause the setting of DIVIDE_BY_ZERO.
Note: It is intended that implementations be free to use code motion
techniques.
The processor need not support INVALID, UNDERFLOW, and INEXACT. If an
exception is not supported, its flag is always clear. The function
IEEE_SUPPORT_FLAG may be used to inquire whether a particular flag is
supported.
The constant USUAL is associated with the list OVERFLOW,
DIVIDE_BY_ZERO, INVALID. The constant ALL is associated with the list
of all the flags.
<<15.3 The rounding modes>>
IEEE 754-1985 specifies four rounding modes:
NEAREST rounds the exact result to the nearest representable value.
TO_ZERO rounds the exact result towards zero to the next
representable value.
UP rounds the exact result towards +infinity to the next
representable value.
DOWN rounds the exact result towards -infinity to the next
representable value.
The function IEEE_GET_ROUNDING_MODE may be used to inquire which
rounding mode is in operation. Its value is one of the above four
or OTHER if the rounding mode does not conform to IEEE 754-1985.
If the processor supports the alteration of the rounding mode during
execution, the subroutine IEEE_SET_ROUNDING_MODE may be used to alter
it. The function IEEE_SUPPORT_ROUNDING may be used to inquire whether
this facility is available for a particular mode.
<<15.4 Halting>>
If the processor supports the ability to control during program
execution whether to abort or continue execution after an exception,
such control may be exercised by invocation of the subroutine
IEEE_HALT. Halting is not precise and may occur any time after the
exception has occurred. The function IEEE_SUPPORT_HALTING may be used
to inquire whether this facility is available.
<<15.5 The floating point status>>
The values of all the supported flags for exceptions, rounding mode,
and halting may be saved in a scalar variable of type
TYPE(IEEE_FSR_VALUE) with the function IEEE_GET_FSR and restored with
the subroutine IEEE_SET_FSR. There are no facilities for finding the
values of particular flags held within such a variable.
Note. Some processors hold all these flags in a floating point
status register that can be saved and restored as a whole much
faster than all individual flags can be saved and restored. These
procedures are provided to exploit this feature.
<<15.6 Exceptional values>>
IEEE 754-1985 specifies the following exceptional floating point values:
Denormalized values have very small absolute values and lowered precision.
Infinite values are created by overflow or division by zero.
Not-a-Number (NaN) values are undefined values or values created by
an invalid operation.
The functions IEEE_IS_DENORMAL, IEEE_IS_INF, IEEE_IS_NAN are provided
to test whether a value is denormalized, infinite, or Nan. The
functions IEEE_INFINITY and IEEE_NAN are provided to generate an
infinity or a Nan. The function IEEE_IS_NORMAL is provided to test
whether a value is none of these. The functions IEEE_SUPPORT_DENORMAL,
IEEE_SUPPORT_INF, and IEEE_SUPPORT_NAN may be used to inquire whether
this facility is available for a particular kind of real.
<<15.7 IEEE arithmetic>>
The function IEEE_DATATYPE may be used to inquire whether IEEE
arithmetic is available for a particular kind of real. Complete
conformance with the IEEE standard is not required, but the normalized
numbers must be exactly those of IEEE single, IEEE double, or the
corresponding quad with an exponent width of 15 bits and a mantissa
width of 112 bits; the arithmetic operators must be implemented with at
least one of the IEEE rounding modes; and the functions copysign,
scalb, logb, nextafter, and unordered must be provided by the functions
IEEE_COPYSIGN, IEEE_SCALB, IEEE_LOGB, IEEE_NEXTAFTER, and
IEEE_UNORDERED.
IEEE 754-1985 specifies a square root function that returns -0 for the
square root of -0. The function IEEE_SQRT is provided for this and the
function IEEE_SUPPORT_SQRT may be used to inquire whether this facility
is available for a particular kind of real.
The inquiry function SUPPORT_ALL is provided to inquire whether the
processor supports all the IEEE facilities defined in this standard for
all kinds of reals.
<<15.8 Tables of the procedures>>
In this section, the procedures are tabulated with the names of their
arguments and a short description.
<<15.8.1 Inquiry functions>>
IEEE_DATATYPE (X) Inquire if the processor supports IEEE arithmetic
for reals of the same kind type parameter as the argument.
IEEE_SUPPORT_ALL() Inquire if processor supports all the IEEE facilities
defined in this standard for all kinds of reals.
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_FLAG(FLAG) Inquire if the processor supports an exception
or set of exceptions.
IEEE_SUPPORT_HALTING() Inquire if the processor supports the ability
to control during program execution whether to abort or continue
execution after an exception.
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_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_ROUNDING (RND_VALUE) Inquire if processor supports
changing to a particular rounding mode for IEEE kinds of reals.
IEEE_SUPPORT_SQRT(X) Inquire if the processor supports IEEE square root
for reals of the same kind type parameter as the argument.
<<15.8.2 Scalar functions>>
IEEE_GET_FLAG(FLAG) Get an exception flag.
IEEE_GET_FSR() Save the current values of the set of flags that define
the current floating point status, including all the exception
flags.
IEEE_GET_ROUNDING_MODE() Save the current IEEE rounding mode.
<<15.8.3 Elemental functions>>
IEEE_COPYSIGN(X,Y) IEEE copysign function.
IEEE_INFINITY(X) Generate IEEE infinity with the same sign as X.
IEEE_IS_DENORMAL(X) Determine if value is IEEE denormalized.
IEEE_IS_INF(X) Determine if value is IEEE Infinity.
IEEE_IS_NAN(X) Determine if value is IEEE Not-a-Number.
IEEE_IS_NORMAL(X) Whether a value is normal, that is, neither an
Infinity, a NaN, nor denormalized.
IEEE_LOGB(X) Unbiased exponent in the IEEE floating point format.
IEEE_NAN(X) Generate IEEE Not-a-Number.
IEEE_NEXTAFTER(X,Y) Returns the next representable neighbor of X in the
direction toward Y.
IEEE_SCALB (X,I) Returns X * 2**I.
IEEE_SQRT(X) Square root as defined in the IEEE standard.
IEEE_UNORDERED(X,Y) IEEE unordered function. True if X or Y is a NaN
and false otherwise.
<<15.8.4 Subroutines>>
IEEE_CLEAR_FLAGS (FLAG1 [, FLAG2, ...]) Clear one or more exception flags.
IEEE_HALT(FLAG,HALTING) Controls continuation or halting on exceptions.
IEEE_SET_FLAGS(FLAG1 [, FLAG2, ...]) Set one or more exception flags.
IEEE_SET_FSR(FSR_VALUE) Restore the values of the set of flags that
define the the floating point status.
IEEE_SET_ROUNDING_MODE(RND_VALUE) Set the current IEEE rounding mode.
<<15.9 Specifications of the procedures>>
IEEE_CLEAR_FLAGS (FLAG1 [, FLAG2, ...])
Description. Clear one or more exception flags.
Class. Subroutine.
Arguments. Each argument shall be of type TYPE(IEEE_FLAG) and with
value OVERFLOW, DIVIDE_BY_ZERO, INVALID, USUAL, UNDERFLOW,
INEXACT, or ALL, which are named PARAMETERS
defined in the module. They are INTENT(IN) arguments and
specify the IEEE flags to be cleared. Specifying USUAL is
equivalent to specifying OVERFLOW, DIVIDE_BY_ZERO,
INVALID. Specifying ALL is equivalent to specifying all
the flags.
Example. CALL IEEE_CLEAR_FLAGS(OVERFLOW,UNDERFLOW) clears the
overflow and underflow flags.
IEEE_COPYSIGN(X,Y)
Description. IEEE copysign function.
Class. Elemental function.
Arguments. The arguments shall be of type real and such that both
IEEE_DATATYPE(X) and IEEE_DATATYPE(Y) have the value true.
Result Characteristics. Same as X.
Result Value. The result has the value of X with the sign of Y.
This is true even for IEEE special values, such as Nan and Inf (on
processors supporting such values).
Examples. The value of COPYSIGN(X,1.0) is ABS(X) even when X is
NaN. The value of COPYSIGN(-X,1.0) is X copied with its sign
reversed, not 0-x; the distinction is germane when X is +0, -0, or
NaN.
IEEE_DATATYPE(X)
Description. Inquire if the processor supports IEEE arithmetic
for reals of the same kind type parameter as the argument.
Class. Inquiry function.
Argument. X shall be scalar and of type real.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor supports
IEEE arithmetic for real and complex variables of the same kind
type parameter as X; otherwise, it has the value false. If the
result is true, the normalized numbers shall be exactly those of
IEEE single, IEEE double, or the corresponding quad with an
exponent width of 15 bits and a mantissa width of 112 bits;
the arithmetic operators shall be implemented with at least
one of the IEEE rounding modes; and the functions IEEE_IS_NORMAL,
IEEE_COPYSIGN, IEEE_SCALB, IEEE_LOGB, IEEE_NEXTAFTER, and
IEEE_UNORDERED shall be provided.
Example. IEEE_DATATYPE(1.0) has the value .TRUE. if default reals
are implemented as in the IEEE standard except that underflowed
values flush to 0.0 instead of being denormal.
IEEE_GET_FLAG(FLAG)
Description. Get an exception flag.
Class. Transformational function.
Argument. FLAG shall be scalar and of type TYPE(IEEE_FLAG). Its value
shall be OVERFLOW, DIVIDE_BY_ZERO, INVALID, USUAL, UNDERFLOW,
INEXACT, or ALL, which are named PARAMETERS defined in the
module. It specifies the IEEE flag to be inspected.
Result Characteristics. Default logical scalar.
Result Value.
Case (i): If the value of FLAG is not USUAL or ALL, the result value
is true if the exception flag is set and is false otherwise.
Case (ii): If the value of FLAG is USUAL, the result value is true if
INVALID, OVERFLOW, or DIVIDE_BY_ZERO is set and is
false otherwise.
Case (iii): If the value of FLAG is ALL, the result value is true if
any flag is set and is false otherwise.
Example. The value of IEEE_GET_FLAG(OVERFLOW) is true if the
overflow flag is set and is false if it is clear.
IEEE_GET_FSR()
Description. Save the current values of the set of flags that define
the current floating point status, including all the exception
flags.
Class. Transformational function.
Arguments. None.
Result Characteristics. Scalar of type TYPE(IEEE_FSR_VALUE).
Result Value. The result value is that of the floating point
status.
Example. To store all the exceptions flags, do a calculation
involving exception handling, and restore them later:
USE IEEE_ARITHMETIC
TYPE(IEEE_FSR_VALUE) FSR_VALUE
:
FSR_VALUE = IEEE_GET_FSR() ! Store the flags
CALL IEEE_CLEAR_FLAGS (ALL)
: ! calculation involving exception handling
CALL IEEE_SET_FSR(FSR_VALUE) ! Restore the flags
Note. The result can be used only in an IEEE_SET_FSR invocation.
IEEE_GET_ROUNDING_MODE()
Description. Save the current IEEE rounding mode.
Class. Transformational function.
Arguments. None.
Result Characteristics. Scalar of type TYPE(IEEE_RND_VALUE).
Result Value. The result value is that of the floating point
rounding mode, and may be NEAREST, TO_ZERO, UP, DOWN,
or OTHER.
Example. To store the rounding mode, do a calculation
with round to nearest, and restore the rounding mode later:
USE IEEE_ARITHMETIC
TYPE(IEEE_RND_VALUE) RND_VALUE
:
RND_VALUE = IEEE_GET_ROUNDING_MODE() ! Store the rounding mode
CALL IEEE_SET_ROUNDING_MODE (NEAREST)
: ! calculation with round to nearest
CALL IEEE_SET_ROUNDING_MODE(RND_VALUE) ! Restore the rounding mode
Note. The result can legally be used only in an IEEE_SET_ROUNDING_MODE
invocation.
IEEE_HALT(FLAG,HALTING)
Description. Controls continuation or halting on exceptions.
Class. Subroutine. IEEE_SUPPORT_HALTING() shall be true.
Arguments.
FLAG shall be scalar and of type TYPE(IEEE_FLAGS). It is of
INTENT(IN) and shall have one of the values: OVERFLOW,
DIVIDE_BY_ZERO, INVALID, USUAL, UNDERFLOW, INEXACT, or ALL,
which are named PARAMETERS defined in the module. USUAL
specifies the exceptions INEXACT, OVERFLOW, and
DIVIDE_BY_ZERO. ALL specifies all exceptions.
HALTING shall be scalar and of type logical. It is of
INTENT(IN). If the value is true, the exception or
exceptions specified by FLAG will cause halting.
Otherwise, execution will continue after these
exceptions.
Example. CALL IEEE_HALT(DIVIDE_BY_ZERO,.TRUE.) causes halting after
a divide_by_zero exception.
Note: Halting is not precise and may occur some time after the
exception has occurred.
IEEE_INFINITY(X)
Description. Generate IEEE infinity.
Class. Elemental function.
Argument. X shall be of type real and such that IEEE_SUPPORT_INF(X)
has the value true.
Result Characteristics. Same as X.
Result Value. The result value is the IEEE infinity with the sign of X.
Example. IEEE_INFINITY(-1.0) has the value -Infinity.
IEEE_IS_DENORMAL(X)
Description. Whether a value is IEEE denormalized.
Class. Elemental function.
Argument. X shall be of type real and such that
IEEE_SUPPORT_DENORMAL(X) has the value true.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the value of X is
IEEE denormalized; otherwise, it has the value false.
Example. IEEE_IS_DENORMAL(TINY(X)/2) has the value true.
IEEE_IS_INF(X)
Description. Whether a value is IEEE Infinity.
Class. Elemental function.
Argument. X shall be of type real and such that
IEEE_SUPPORT_INF(X) has the value true.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the value of X is an
IEEE Infinity (either positive or negative); otherwise, it has the
value false.
Example. IEEE_IS_INF(IEEE_INFINITY(-1.0)) has the value true.
Note: the sign of X can be determined directly via comparison (e.g.
(X > 0) is true for +Inf) or indirectly via SIGN or COPYSIGN.
IEEE_IS_NAN(X)
Description. Whether a value is IEEE Not-a-Number.
Class. Elemental function.
Argument. X shall be of type real and such that
IEEE_SUPPORT_NAN(X) has the value true.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the value of X is an
IEEE NaN; otherwise, it has the value false.
Example. IEEE_IS_NAN(IEEE_SQRT(-1.0)) has the value true.
IEEE_IS_NORMAL(X)
Description. Whether a value is normal, that is, neither an
Infinity, a NaN, nor denormalized.
Class. Elemental function.
Argument. X shall be of type real and such that IEEE_DATATYPE(X) has
the value true.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the value of X is
neither an IEEE Infinity (either positive or negative), an IEEE
NaN, nor denormalized; otherwise, it has the value false.
Example. IEEE_IS_NORMAL(IEEE_SQRT(-1.0)) has the value false.
Note: While one could code this via a combination of the primitives,
it is often significantly more efficient to just test it once. It also
makes for more readable code.
IEEE_LOGB(X)
Description. Unbiased exponent in the IEEE floating point format.
Class. Elemental function.
Argument. X shall be of type real and such that IEEE_DATATYPE (X)
has the value true.
Result Characteristics. Default integer.
Result Value.
Case (i): If the value of X is neither zero, infinity, nor NaN,
the result has the value of the unbiased exponent of X.
Case(ii): If X==0, the result is -HUGE(X).
Case(iii): If X has an infinite value, the result is +Inf.
Case(iv): If X has a NaN value, the result is X.
Note: Cases (iii) and (iv) cannot occur on processors that lack Inf and
NaN values.
Example. IEEE_LOGB(-1.1) has the value 0.
Comment to Keith: I have made the result integer.
IEEE_NAN(X)
Description. Generate IEEE Not-a-Number.
Class. Elemental function.
Argument. X shall be of type real and such that IEEE_SUPPORT_NAN(X)
has the value true.
Result Characteristics. Same as X.
Result Value. The result value is an IEEE signaling NaN.
Example. IEEE_NAN(X) is an IEEE signaling NaN.
Note: IEEE defines an entire family of values as NaN. A processor may
choose to allow users to specify the bit pattern. If so, the value of
X is appropriately transformed to be a signaling NaN with all the low
order bits taken from the value of X.
Comment: one can argue that IEEE_NAN and IEEE_INF are not needed,
because users can code expressions to generate them. However,
experience shows that optimizers often make it difficult to do this
easily and portably. It is better to have these functions in the
Standard to promote portability, without impacting optimization.
IEEE_NEXTAFTER(X,Y)
Description. Returns the next representable neighbor of X in the
direction toward Y.
Class. Elemental function.
Arguments. The arguments shall be of type real and such that
IEEE_DATATYPE(X) and IEEE_DATATYPE(Y) have the value true.
Result Characteristics. Same as X.
Result Value.
Case (i): If X == Y, the result is X without any exception ever
being signaled.
Case (ii): If X /= Y, the result has the value of the next
representable neighbor of X in the direction of Y. If
either X or Y is a NaN, the result is one or the
other of the input NaNs. Overflow is signaled when X is
finite but NEXTAFTER(X,Y) is infinite; underflow is
signaled when NEXTAFTER(X,Y) is denormalized; in both
cases, INEXACT is signaled.
Example. The value of IEEE_NEXTAFTER(1.0,2.0) is 1.0+EPSILON(X).
IEEE_SCALB(X,I)
Description. Returns X * 2**I.
Class. Elemental function.
Arguments.
X shall be of type real and such that IEEE_DATATYPE(X)
has the value true.
I shall be of type integer.
Result Characteristics. Same as X.
Result Value.
Case (i): If X * 2**I is within range, the result has this value.
Case (ii): If X * 2**I is too large and IEEE infinities are supported,
the result value is infinity with the sign of X.
Case (iii): If X * 2**I is too large and IEEE infinities are not supported,
the result value is SIGN(HUGE,X).
Example. The value of IEEE_SCALB(1.0,2) is 4.0.
IEEE_SET_FLAGS(FLAG1 [, FLAG2, ...])
Description. Set one or more supported exception flags.
Class. Subroutine.
Arguments. The arguments shall be of type TYPE(IEEE_FLAG) and with
value OVERFLOW, DIVIDE_BY_ZERO, INVALID, USUAL, UNDERFLOW,
INEXACT, or ALL, which are named constants defined in the
module. They are INTENT(IN) arguments and
specify the exception flags to be set. USUAL is equivalent
to the list OVERFLOW, DIVIDE_BY_ZERO, INVALID. ALL is
equivalent to the list of all flags. If a flag is not
supported on the processor, it remains clear.
Example. CALL IEEE_SET_FLAGS(OVERFLOW) sets the overflow flag.
IEEE_SET_FSR(FSR_VALUE)
Description. Restore the values of the set of flags that define the
the floating point status.
Class. Subroutine.
Argument. FSR_VALUE shall be scalar and of type TYPE(IEEE_FSR_VALUE).
Its value shall have been set in a previous invocation of IEEE_GET_FSR.
Example. To store all the exceptions flags, do a calculation
involving exception handling, and restore them later:
USE IEEE_ARITHMETIC
TYPE(IEEE_FSR_VALUE) FSR_VALUE
:
FSR_VALUE = IEEE_GET_FSR() ! Store the flags
CALL IEEE_CLEAR_FLAGS (ALL)
: ! calculation involving exception handling
CALL IEEE_SET_FSR(FSR_VALUE) ! Restore the flags
Note: getting and setting may be expensive operations. It is the
programmer's responsibility to do it when necessary to assure correct
results.
IEEE_SET_ROUNDING_MODE (RND_VALUE)
Description. Set the current IEEE rounding mode.
Class. Subroutine. IEEE_SUPPORT_ROUNDING(RND_VALUE) shall be true.
Argument. RND_VALUE shall be scalar and of type TYPE(IEEE_RND_VALUE).
Its value shall be one of NEAREST, TO_ZERO, UP, and DOWN, which are
named constants of the module.
Example. To store the rounding mode, do a calculation
with round to nearest, and restore the rounding mode later:
USE IEEE_ARITHMETIC
TYPE(IEEE_RND_VALUE) RND_VALUE
:
RND_VALUE = IEEE_GET_ROUNDING_MODE() ! Store the rounding mode
CALL IEEE_SET_ROUNDING_MODE (NEAREST)
: ! calculation with round to nearest
CALL IEEE_SET_ROUNDING_MODE(RND_VALUE) ! Restore the rounding mode
IEEE_SQRT(X)
Description. Square root as defined in the IEEE standard.
Class. Elemental function.
Argument. X shall be of type real and such that IEEE_SUPPORT_SQRT(X)
has the value true.
Result characteristics. Same as X.
Result value. The result has the value of the square-root of X as
defined in the IEEE Standard.
Example. IEEE_SQRT(-0.) has the value -0.
Note: On most IEEE processors, SQRT is probably the same as IEEE_SQRT.
However, processors without adequate hardware support may choose to
return different values in hard cases, in which case a careful
programmer may wish to use IEEE_SQRT for safety.
IEEE_SUPPORT_ALL()
Description. Inquire if processor supports all the IEEE facilities
defined in this standard for all kinds of reals.
Class. Inquiry function.
Arguments. None.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the results of all
the functions IEEE_DATATYPE, IEEE_SUPPORT_DENORMAL, IEEE_SUPPORT_FLAG,
IEEE_SUPPORT_HALTING, IEEE_SUPPORT_INF, IEEE_SUPPORT_NAN,
IEEE_SUPPORT_ROUNDING, and IEEE_SUPPORT_SQRT are always true;
otherwise, it has the value false.
Example. IEEE_SUPPORT_ALL() has the value false if the processor
supports both IEEE and non-IEEE kinds of reals.
IEEE_SUPPORT_DENORMAL(X)
Description. Inquire if the processor supports IEEE gradual underflow
for reals of the same kind type parameter as the argument.
Class. Inquiry function.
Argument. X shall of type real and such that IEEE_DATATYPE(X) has
the value true. It may be scalar or array valued.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor supports
IEEE gradual underflow for real variables of the same kind type
parameter as X; otherwise, it has the value false.
Example. IEEE_SUPPORT_DENORMAL(X) has the value true if the processor
supports gradual underflow for X.
IEEE_SUPPORT_FLAG(FLAG)
Description. Inquire if the processor supports an exception or set
of exceptions.
Class. Inquiry function.
Argument. FLAG shall be scalar and of type TYPE(IEEE_FLAGS). It
shall have one of the values: INEXACT, UNDERFLOW, OVERFLOW, INVALID,
DIVIDE_BY_ZERO, USUAL and ALL. USUAL specifies the exceptions
INEXACT, OVERFLOW, and DIVIDE_BY_ZERO. ALL specifies the exceptions
INEXACT, UNDERFLOW, OVERFLOW, INVALID, and DIVIDE_BY_ZERO.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor supports
detection of the specified exception or exceptions; otherwise, it
has the value false.
Example. IEEE_SUPPORT_FLAG(ALL) has the value true if the processor
supports all the exceptions.
IEEE_SUPPORT_HALTING()
Description. Inquire if the processor supports the ability to control
during program execution whether to abort or continue execution
after an exception.
Class. Inquiry function.
Arguments. None
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor
supports the ability to control during program execution whether to
abort or continue execution an exception; otherwise, it has
the value false.
Example. IEEE_SUPPORT_HALTING() has the value true if the processor
supports halting after an exception.
IEEE_SUPPORT_INF(X)
Description. Inquire if processor supports the IEEE infinity facility
for reals of the same kind type parameter as the argument.
Class. Inquiry function.
Argument. X shall of type real and such that IEEE_DATATYPE(X) has
the value true. It may be scalar or array valued.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor supports
IEEE infinities (positive and negative) for real variables of the
same kind type parameter as X; otherwise, it has the value false.
Example. IEEE_SUPPORT_INF(X) has the value true if the processor
supports IEEE infinities for X.
IEEE_SUPPORT_NAN(X)
Description. Inquire if processor supports the IEEE Not-A-Number
facility for reals of the same kind type parameter as the argument.
Class. Inquiry function.
Argument. X shall of type real and such that IEEE_DATATYPE(X) has
the value true. It may be scalar or array valued.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor
supports IEEE NaNs for real variables of the same kind type
parameter as X; otherwise, it has the value false.
Example. IEEE_SUPPORT_NAN(X) has the value true if the processor
supports IEEE NaNs for X.
IEEE_SUPPORT_ROUNDING(RND_VALUE)
Description. Inquire if processor supports changing to a particular
rounding mode for IEEE kinds of reals.
Class. Inquiry function.
Argument. RND_VALUE shall be of type TYPE(IEEE_RND_VALUE) and value
one of NEAREST, TO_ZERO, UP, and DOWN, which are named constants of
the module.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor
supports the IEEE rounding mode defined by RND_VALUE for IEEE kinds
of reals; otherwise, it has the value false.
Example. IEEE_SUPPORT_ROUNDING(TO_ZERO) has the value true if the
processor supports rounding to zero for IEEE kinds of reals.
IEEE_SUPPORT_SQRT(X)
Description. Inquire if the processor supports IEEE square root
for reals of the same kind type parameter as the argument.
Class. Inquiry function.
Argument. X shall of type real and such that IEEE_DATATYPE(X) has
the value true. It may be scalar or array valued.
Result Characteristics. Default logical scalar.
Result Value. The result has the value true if the processor supports
IEEE square root for real variables of the same kind type parameter
as X; otherwise, it has the value false.
Example. IEEE_SUPPORT_SQRT(X) has the value true if the processor
supports IEEE square root for X.
IEEE_UNORDERED(X,Y)
Description. IEEE unordered function. True if X or Y is a NaN.
and false otherwise.
Class. Elemental function.
Arguments. The arguments shall be of type real and such that
IEEE_DATATYPE(X) and IEEE_DATATYPE(Y) have the value true.
Result Characteristics. Same as X.
Result Value. The result has the value true if X or Y is a NaN
or both are NaNs; otherwise, it has the value false.
Example. IEEE_UNORDERED(0.,IEEE_SQRT(-1.0)) has the value true.
<<15.10 Examples>>
Example 1:
MODULE DOT
! Module for dot product of two real arrays of rank 1.
USE IEEE_ARITHMETIC
LOGICAL MATRIX_ERROR
INTERFACE OPERATOR(.dot.)
MODULE PROCEDURE MULT
END INTERFACE
CONTAINS
REAL FUNCTION MULT(A,B)
REAL, INTENT(IN) :: A(:),B(:)
INTEGER I
LOGICAL OLD_OVERFLOW
IF (SIZE(A)/=SIZE(B)) THEN
MATRIX_ERROR = .TRUE.
RETURN
END IF
OLD_OVERFLOW = IEEE_GET_FLAG(OVERFLOW)
CALL IEEE_CLEAR_FLAGS(OVERFLOW)
MULT = 0.
DO I = 1, SIZE(A)
MULT = MULT + A(I)*B(I)
END DO
IF (IEEE_GET_FLAG(OVERFLOW) ) THEN
MATRIX_ERROR = .TRUE.
ELSE
IF(OLD_OVERFLOW) CALL IEEE_SET_FLAGS(OVERFLOW)
END IF
END FUNCTION MULT
END MODULE DOT
This module provides the dot product of two real arrays of rank 1. If
the sizes of the arrays are different, an immediate return occurs with
MATRIX_ERROR true. If OVERFLOW occurs during the actual calculation,
the overflow flag will be set and MATRIX_ERROR will be true.
Example 2:
USE IEEE_ARITHMETIC
TYPE(IEEE_FSR_VALUE) FSR_VALUE
FSR_VALUE = IEEE_GET_FSR()
CALL IEEE_CLEAR_FLAGS (OVERFLOW, INVALID, DIVIDE_BY_ZERO)
! First try the "fast" algorithm for inverting a matrix:
MATRIX1 = FAST_INV (MATRIX) ! This must not alter MATRIX.
IF (IEEE_GET_FLAG(OVERFLOW) .OR. IEEE_GET_FLAG(INVALID) .OR. &
IEEE_GET_FLAG(DIVIDE_BY_ZERO) ) THEN
! "Fast" algorithm failed; try "slow" one:
CALL IEEE_CLEAR_FLAGS (OVERFLOW, INVALID, DIVIDE_BY_ZERO)
MATRIX1 = SLOW_INV (MATRIX)
IF (IEEE_GET_FLAG(OVERFLOW) .OR. IEEE_GET_FLAG(INVALID) .OR. &
IEEE_GET_FLAG(DIVIDE_BY_ZERO) ) THEN
WRITE (*, *) 'Cannot invert matrix'
STOP
END IF
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 3:
USE IEEE_ARITHMETIC
:
! First try a fast algorithm for inverting a matrix.
CALL IEEE_CLEAR_FLAGS (OVERFLOW)
DO K = 1, N
:
IF (IEEE_GET_FLAG(OVERFLOW)) EXIT
END DO
IF (IEEE_GET_FLAG(OVERFLOW)) THEN
! Alternative code which knows that K-1 steps have executed normally.
:
END IF
Here the code for matrix inversion is in line and the transfer is made more
precise by adding extra tests of the flag.
Example 4:
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_GET_FLAG(OVERFLOW)
OLD_UNDERFLOW = IEEE_GET_FLAG(UNDERFLOW)
CALL IEEE_CLEAR_FLAGS(UNDERFLOW, OVERFLOW)
! Try a fast algorithm first
HYPOT = SQRT( X**2 + Y**2 )
IF (IEEE_GET_FLAG(UNDERFLOW) .OR. IEEE_GET_FLAG(OVERFLOW) ) THEN
CALL IEEE_CLEAR_FLAGS(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_SET_FLAGS(OVERFLOW)
IF(OLD_UNDERFLOW) CALL IEEE_SET_FLAGS(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.