| Fortran 95 Interval Arithmetic Programming Reference |    
|
| Fortran 95 Interval Arithmetic Programming Reference |
|
| C H A P T E R 1 |
|
Using Interval Arithmetic With f95 |
Interval arithmetic is a system for computing with intervals of numbers. Because interval arithmetic always produces intervals that contain the set of all possible result values, interval algorithms have been developed to perform surprisingly difficult computations. For more information on interval applications, see the Interval Arithmetic Readme.
Since the inception of interval arithmetic, interval algorithms that produce narrow-width results have been developed, and the syntax and semantics for interval language support have been designed. However, relatively little progress has been made in providing commercially available and supported interval compilers. With one exception (M77 Minnesota FORTRAN 1977 Standards Version Edition 1), interval systems have been based on pre-processors, C++ classes, or Fortran 90 modules. The goals of intrinsic compiler support for interval data types in f95 are:
Sun Studio Fortran 95 interval support is a significant extension to Fortran.
The goal of intrinsic INTERVAL support in f95 is to stimulate development of commercial interval solver libraries and applications by providing program developers with:
Support and features are components of implementation quality. Not all possible quality of implementation features have been implemented. Throughout this book, various unimplemented quality of implementation opportunities are described. Additional suggestions from users are welcome.
As a consequence of evaluating any interval expression, a valid interval-supporting compiler must produce an interval that contains the set of all possible results. The set of all possible results is called the containment set (cset) of the given expression. The requirement to enclose an expression's cset is the containment constraint of interval arithmetic. The failure to satisfy the containment constraint is a containment failure. A silent containment failure (with no warning or documentation) is a fatal error in any interval computing system. By satisfying this single constraint, intervals provide otherwise unprecedented computing quality.
Given the containment constraint is satisfied, implementation quality is determined by the location of a point in the two-dimensional plane whose axes are runtime and interval width. On both axes, small is better. How to trade runtime for interval width depends on the application. Both runtime and interval width are obvious measures of interval-system quality. Because interval width and runtime are always available, measuring the accuracy of both interval algorithms and implementation systems is no more difficult than measuring their speed.
The Sun Studio tools for performance profiling can be used to tune interval programs. However, in f95, no interval-specific tools exist to help isolate where an algorithm may gain unnecessary interval width. As described in Section 1.4, Code Development Tools, some interval dbx and global program checking (GPC) support are provided. Adding additional interval-specific code development and debugging tools are quality of implementation opportunities.
All the normal language and compiler quality of implementation opportunities exist for intervals, including rapid execution and ease-of-use.
Valid interval implementation systems include a new additional quality of implementation opportunity: Minimize the width of computed intervals while always satisfying the containment constraint.
If an interval's width is as narrow as possible, it is said to be sharp. For a given floating-point precision, an interval result is sharp if its width is as narrow as possible.
The following can be said about the width of intervals produced by the f95 compiler:
By providing compiler optimization and hardware instruction support, INTERVAL operations are not necessarily slower than their REAL floating-point counterparts. In f95, the following can be said about the speed of intrinsic interval operators and mathematical functions:
The intrinsic INTERVAL data type in Fortran facilitates interval code development, testing, and execution. To make interval code transparent (easy to write and read), interval syntax and semantics have been added to Fortran. User acceptance will ultimately determine which interval features are added to standard Fortran.
By introducing intervals as an intrinsic data type to Fortran, all of the applicable syntax and semantics of the Fortran language become immediately available. Sun Studio Fortran 95 includes the following interval-specific Fortran extensions:
For examples and more information on these and other intrinsic interval functions, see CODE EXAMPLE 1-11 through CODE EXAMPLE 1-14 and Section 2.10.4.5, Intrinsic Functions.
Chapter 2 contains detailed descriptions of these and other interval features.
The examples in this section are designed to help new interval programmers to understand the basics and to quickly begin writing useful interval code. Modifying and experimenting with the examples is strongly recommended.
All code examples in this book are contained in the directory:
/opt/SUNWspro/examples/intervalmath/docExamples
The name of each file is cen-m.f95, where n is the chapter in which the example occurs, and m is the number of the example. Additional interval examples are contained in the directory:
/opt/SUNWspro/examples/intervalmath/general
The following f95 command-line macro is the simplest way to invoke recognition of INTERVAL data types as intrinsic and to control INTERVAL expression processing:
For intrinsic INTERVAL data types to be recognized by the compiler, either -xia or -xinterval must be entered in the f95 command line.
All command-line options that interact with intervals are described in Section 2.3.3, Interval Command-Line Options. Widest-need and strict expression processing are described in Section 2.3, INTERVAL Arithmetic Expressions.
The simplest command-line invocation of f95 with interval support is shown in CODE EXAMPLE 1-1.
Unless explicitly stated otherwise, all code examples are compiled using the -xia command-line option. The -xia or -xinterval command-line option is required to use the interval extensions to f95.
CODE EXAMPLE 1-1 is the interval equivalent of "hello world."
math% cat ce1-1.f95 PRINT *, "[2, 3] + [4, 5] = ", [2, 3] + [4, 5] ! line 1 END math% f95 -xia ce1-1.f95 math% a.out [2, 3] + [4, 5] = [6.0,8.0] |
CODE EXAMPLE 1-1 uses list-directed output to print the labeled sum of the intervals [2, 3] and [4, 5].
The INTERVAL declaration statement performs the same functions for INTERVAL data items as the REAL, INTEGER, and COMPLEX declarations do for their respective data items. The default INTERVAL kind type parameter value (KTPV) is twice the default INTEGER KTPV. This permits any default INTEGER to be exactly represented using a degenerate default INTERVAL. See Section 1.3.7, Default Kind Type Parameter Value (KTPV) for more information.
CODE EXAMPLE 1-2 uses INTERVAL variables and initialization to perform the same operation as CODE EXAMPLE 1-1.
In line 1, the variables, X and Y are declared to be default type INTERVAL variables and are initialized to [2, 3] and [4, 5], respectively. Line 2 uses list-directed output to print the labeled interval sum of X and Y.
Full support for reading and writing intervals is provided. Reading and writing INTERVAL and COMPLEX data items are similar. Intervals use square brackets, instead of parentheses as delimiters.
In f95 the input conversion process constructs a sharp interval that contains the input decimal value. If the value is machine representable, the internal machine approximation is degenerate. If the value is not machine representable, an interval having width of 1-ulp (unit-in-the-last-place of the mantissa, pronounced "ulp") is constructed.
The simplest way to read and print INTERVAL data items is with list-directed input and output.
CODE EXAMPLE 1-3 is a simple tool to help users become familiar with interval arithmetic and single-number INTERVAL input/output using list-directed READ and PRINT statements. Complete support for formatted INTERVAL input/output is provided, as described in Section 2.10.2, Input and Output.
|
Note - The interval containment constraint requires that directed rounding be used during input and output. With single-number input followed immediately by single-number output, a decimal digit of accuracy can appear to be lost. In fact, the width of the input interval is increased by at most 1-ulp, when the input value is not machine representable. See Section 1.3.5, Single-Number Input/Output and CODE EXAMPLE 1-6 |
|
Note - The empty interval is supported in f95. The empty interval can be entered as "[empty]". Infinite interval endpoints are also supported, as described in Section 2.10.2.1, External Representations and illustrated in CODE EXAMPLE 2-37. |
One of the most frustrating aspects of reading interval output is comparing interval infima and suprema to count the number of digits that agree. For example, CODE EXAMPLE 1-4 and CODE EXAMPLE 1-5 shows the interval output of a program that generates different random width INTERVAL data.
|
Note - Only program output is shown in CODE EXAMPLE 1-4 and CODE EXAMPLE 1-5. The code that generates the output is included with the examples locatedat http://developer.sun.com/prodtech/cc/reference/codesamples/ |
Compare the output readability in CODE EXAMPLE 1-4 with CODE EXAMPLE 1-5.
Because reading and interactively entering interval data can be tedious, a single-number interval format is introduced. The single-number convention is that any number not contained in brackets is interpreted as an interval whose lower and upper bounds are constructed by subtracting and adding 1 unit to the last displayed digit.
Thus during interval input and output,
[2.34499, 2.34501] = 2.34500 + [-1, +1]uld
where [-1, +1]uld means that the interval [-1, +1] is added to the last digit of the preceding number. The subscript, uld, is a mnemonic for "unit in the last digit."
|
Note - The single number input/output representation is not used to represent INTERVAL literal constants in f95 code. |
To represent a degenerate interval, a single number can be enclosed in square brackets. For example,
[2.345] = [2.345, 2.345] = 2.345000000000.....
This convention is used both for single-number input/output and to represent degenerate literal INTERVAL constants in Fortran code. Thus, type [0.1] to enter an exact decimal number, even though 0.1 is not machine representable.
During input to a program, both [0.1,0.1] and [0.1] represents the point, 0.1. However, the single-number input/output value 0.1 represents the interval
In the single-number display format, trailing zeros are significant. See Section 2.10.2, Input and Output for more information.
Intervals can always be entered and displayed using the traditional [inf, sup] display format. In addition, a single number in square brackets denotes a point. For example, on input, [0.1] is interpreted as the number 1/10. To guarantee containment, directed rounding is used to construct an internal approximation that is known to contain the number 1/10.
CODE EXAMPLE 1-6 notes:
The f95 compiler contains the following INTERVAL-specific statements, expressions, and extensions:
In f95 the default INTEGER KTPV is KIND(0) = 4. To represent any default INTEGER with a degenerate default INTERVAL requires the default INTERVAL KTPV, KIND([0]), to be 2*KIND(0) = 8. Choosing 8 for the default INTERVAL KTPV is also done because:
TABLE 1-1 notes:
1. The letters a and b are placeholders for literal decimal constants, such as 0.1 and 0.2.
2. A single decimal constant contained in square brackets denotes a degenerate INTERVAL constant. The same convention is used in input/output.
3. Let expr stand for any Fortran arithmetic expression, whether or not it contains items of type INTERVAL. An assignment statement, V = expr, evaluates the expression, expr, and assigns the resulting value to V. Mixed-mode INTERVAL expressions are not permitted under the -xia=strict command line option. Under the -xia or -xia=widestneed option, mixed-mode expressions are correctly evaluated using widest-need expression processing. Before expression evaluation under widest-need, all integer and floating-point data items are promoted to containing intervals with the largest KTPV found anywhere in the expression, including, V. For details, see Section 2.3.2, Value Assignment.
4. Interval input/output support is designed to provide flexibility, readability, and ease of code development. The most important new edit descriptor is Y, which is used to read and display intervals using the single-number interval format. For a complete description of all edit descriptors that can process intervals, see Section 2.10.2, Input and Output.
The INTERVAL assignment statement assigns the value of an interval expression, denoted by the placeholder expr, to an INTERVAL variable, array element, array, array section, or structure component V. The syntax is:
where V must have an INTERVAL type, and expr denotes any non-COMPLEX numeric expression. Under widest-need expression processing, the expression expr need not be an INTERVAL expression. Under strict expression processing, expr must be an INTERVAL expression with the same KTPV as V.
Gracefully handling mixed-type INTERVAL expressions is an important ease-of-use feature, because it facilitates writing transparent (easy to understand) mathematical expressions.
Mixed-type INTERVAL expressions are supported to make writing and reading interval code no more difficult than it is for REAL code. The interval containment constraint is satisfied in mixed-mode expressions using either widest-need or strict expression processing.
Computing narrow-width interval results is facilitated if the width of INTERVAL constants is dynamically defined by expression context, as described in Section 2.3, INTERVAL Arithmetic Expressions. In mixed-KTPV expressions, shown in CODE EXAMPLE 1-7, dynamically increasing the KTPV of INTERVAL variables can also decrease the width of INTERVAL expression results.
In line 3, KTPVmax = KIND(Z) = 8. This value is used to promote the KTPV of X and Y to 8 before computing their product and storing the result in Z1.
These steps are shown explicitly in the equivalent strict code in line 4.
The process of scanning a statement to determine the maximum KTPV and performing the necessary promotions, is called widest-need expression processing, see Section 2.3, INTERVAL Arithmetic Expressions.
For syntax and semantics of the intrinsic INTERVAL constructor functions, see Section 2.9, Extending Intrinsic INTERVAL Operators.
If the widest-need principle is used with both KTPVs and data types, mixed-mode (type and KTPV) INTERVAL expressions can be safely and predictably evaluated. For example, in CODE EXAMPLE 1-8, the expression for Y1 in line 3 is an interval expression, because X and Y1 are INTERVAL variables.
To guarantee containment, a containing interval must be used in place of a real approximation to the constant 0.1. However, KTPVmax = 16, because KIND(X) = 16. Therefore, the INTERVAL constant [0.1_16], a sharp KTPV = 16 interval containing the exact value, 1/10, is used to update X. Finally, the result is converted to a KTPV = 4 containing interval and assigned to Y1. Line 4 contains the equivalent strict code. Under strict expression processing, neither mixed-type nor mixed-KTPV expressions are permitted.
The logical steps in widest-need expression processing are:
1. Scan the entire statement, including the left-hand side, for any INTERVAL data items.
The presence of any INTERVAL constants, variables, or intrinsic functions, makes the expression's type INTERVAL.
2. Scan the INTERVAL expressions for KTPVmax, based on the KTPV of each INTERVAL, REAL, INTEGER, constant, or variable.
|
Note - Integers are converted to intervals with twice their KTPV so all integer values can be exactly represented. |
3. Promote all variables and constants to intervals with KTPVmax.
5. Convert the result to a lower KTPV if needed to match the left-hand side's KTPV.
6. Assign the resulting value to the left-hand side.
These steps guarantee that mixed-mode INTERVAL expression processing satisfies the containment constraint and efficiently produces reasonably narrow interval results.
Mixed-mode INTERVAL expression evaluation using widest-need expression processing is supported by default with the -xia command-line flag. Using -xia=strict eliminates any automatic type conversions to intervals and any automatic KTPV increases of INTERVAL variables. In strict mode, all interval type and precision conversions must be explicitly coded.
Writing arithmetic expressions that contain INTERVAL data items is simple and straightforward. Except for INTERVAL literal constants and intrinsic INTERVAL-specific functions, INTERVAL expressions look like REAL arithmetic expressions. In particular, with widest-need expression processing, REAL and INTEGER variables and literal constants can be freely used anywhere in an INTERVAL expression, such as in CODE EXAMPLE 1-9.
Because X, the variable to which the assignment is made in line 5, is an INTERVAL, the following steps are taken before evaluating the expression 0.1*A/N:
1. The literal constant 0.1 is converted to the default INTERVAL variable containing the degenerate interval [0.1].
While not required in a valid interval system implementation, Sun Studio Fortran 95 performs sharp data conversions. For example, the internal approximation of [0.1] is 1-ulp wide.
2. The REAL variable A is converted to the degenerate interval [5].
3. The INTEGER variable N is converted to the degenerate interval [3].
The expression [0.1] × [5]/[3] is evaluated using interval arithmetic. The above steps are part of widest-need expression processing, which is required to satisfy the containment constraint when evaluating mixed-mode INTERVAL expressions. See Section 1.3.9, Mixed-Type Expression Evaluation.
An INTERVAL assignment statement must satisfy one requirement: the variable to which the assignment is made must be an INTERVAL variable, array element, array, array section, or structure component. For more information on the widest-need processing mode, see Section 2.3.1, Mixed-Mode INTERVAL Expressions.
Because the interval system implemented in Sun Studio Fortran 95 is closed, if any INTERVAL expression fails to produce a valid interval result, it is a compiler error that should be reported. See Section 1.4, Code Development Tools for information on how to report a suspected error and Section 1.5.1, Known Containment Failures for a list of known errors.
Ordering intervals is more complicated than ordering points. Testing whether 2 is less than 3 is unambiguous. With intervals, while the interval [2,3] is certainly less than the interval [4,5], what should be said about [2,3] and [3,4]?
Three different classes of INTERVAL relational operators are implemented:
For a certainly-relation to be true, every element of the operand intervals must satisfy the relation. A possibly-relation is true if it is satisfied by any elements of the operand intervals. The set-relations treat intervals as sets. The three classes of INTERVAL relational operators converge to the normal relational operators on points if both operand intervals are degenerate.
To distinguish the three operator classes, the normal two-letter Fortran relation mnemonics are prefixed with the letters C, P, or S. In f95 the set operators .SEQ. and .SNE. are the only operators for which the point defaults (.EQ. or == and .NE. or /=) are supported. In all other cases, the relational operator class must be explicitly identified, as for example in:
See Section 2.4, Intrinsic Operators for the syntax and semantics of all INTERVAL operators.
The following program demonstrates the use of a set-equality test.
math% cat ce1-10.f95 INTERVAL :: X = [2, 3], Y = [4, 5] ! Line 1 IF(X+Y .SEQ. [6, 8]) PRINT *, "Check." ! Line 2 END math% f95 -xia ce1-10.f95 math% a.out Check. |
Line 2 uses the set-equality test to verify that X+Y is equal to the interval [6, 8].
IF(X+Y == [6, 8]) PRINT *, "Check." ! line 2
Use CODE EXAMPLE 1-11 and CODE EXAMPLE 1-12 to explore the result of INTERVAL-specific relational operators.
CODE EXAMPLE 1-12 demonstrates the use of the INTERVAL-specific operators ed in TABLE 1-2.
.IH. |
||
.IX. |
||
.DJ. |
||
.IN. |
||
.INT. |
||
.PSB. |
||
.PSP. |
||
.SB. |
||
.SP. |
A variety of intrinsic INTERVAL-specific functions are provided. See Section 2.10.4.5, Intrinsic Functions. Use CODE EXAMPLE 1-13 to explore how intrinsic INTERVAL functions behave.
Every Fortran intrinsic function that accepts REAL arguments has an interval version. See Section 2.10.4.5, Intrinsic Functions. Use CODE EXAMPLE 1-14 to explore how some intrinsic functions behave.
Information on interval code development tools is available online. See the Interval Arithmetic Readme for a list of interval web sites and other online resources.
To report a suspected interval error, send email to
Include the following text in the Subject line of the email message:
FORTEDEV "7.0 mm/dd/yy" Interval
where mm/dd/yy is the month, day, and year.
In Sun Studio, interval data types are supported by dbx to the following extent:
For additional details on dbx functionality, see Debugging a Program With dbx.
Global program checking (GPC) in Sun Studio Fortran 95 detects one interval-specific error: INTERVAL type mismatches in user-supplied routine calls. CODE EXAMPLE 1-15 shows an example of GPC detecting an INTERVAL type mismatch.
The following libraries contain intrinsic INTERVAL routines.
There is limited legacy interval Fortran code with which to contend. Until language syntax and semantics are standardized, different providers of interval compiler support will inevitably diverge. The standardization process will be facilitated if users provide feedback regarding the most favored INTERVAL syntax and semantics. Comments can be sent to the email alias ed in the Interval Arithmetic Readme.
The representation of intervals in binary files will change as compilers supporting narrower interval systems are made available.
In this release, the -autopar compiler option has no effect on loops containing interval arithmetic operations. These loops are not automatically parallelized. The -explicitpar compiler option must be used to parallelize loops marked with explicit parallelization directives.
The following code samples list interval-specific error messages. Each code sample includes the error message and the sample code that produced the error.
Whenever an interval containment failure can occur, a compile-time warning should be issued. There are no know containment failures under widest-need expression processing. In -xia=strict mode, it is possible to violate the containment constraint with an interval ** (integer expression) operation if the integer expression overflows.
Numerical inaccuracies are normally associated with REAL rather than INTEGER expressions. In one respect, INTEGER expressions are harder to detect than REAL expressions. When REAL expressions overflow, an exception is raised and an IEEE infinity is generated. The exception is a warning that overflow has occurred. Infinities tend to propagate in floating-point computations, thereby alerting users of a potential problem. It is also possible to trap on overflow.
When INTEGER expressions overflow, they silently wrap around to some possibly-opposite-signed value. Moreover, the only practical way to detect integer overflow is to perform the inverse operation and test for equality on every integer operation. Integer constant expressions are safe because they are evaluated during compilation where overflow is detected and signalled with a warning message.
Under -xia=widestneed expression processing when the second operand of the ** operator is an integer expression that overflows, the returned interval is guaranteed to contain the correct result. However, the same is not true under -xia=strict processing, because it is not possible to promote integers to intervals prior to evaluating the given expression without widest-need expression processing, The same is true if the second operand of the ** operator is the INTERVAL type conversion routine.
CODE EXAMPLE 1-21 shows that widest-need expression processing is extended to all intrinsic INTEGER operations and functions inside integer expressions in the second operand of the ** operator. This is not true under -xia=strict mode.
This code demonstrates a silent containment failure in -xia=strict mode and the correct interval results in -xia=widestneed mode. For information on the power operator, see Section 2.5, Power Operators X**N and X**Y.
| Fortran 95 Interval Arithmetic Programming Reference | 817-6704-10 |
|
Copyright © 2004, Sun Microsystems, Inc. All rights reserved.
