# Self-Adapting Fortran 77 Machine Constants: Comment on Algorithm 528.

This note discusses user dissatisfaction with the need to uncomment data statements in Algorithm 528, comments on alternative approaches tried by the community, and proposes a solution that is both automatic and safe.Categories and Subject Descriptors: D.3.2 [Programming Languages]: Language Classifications--Fortran 77; G.1.0 [Numerical Analysis]: General--computer arithmetic

General Terms: Algorithms, Languages

Additional Key Words and Phrases: d1mach, machine environment parameters

1. BACKGROUND

For picking stopping tolerances in optimization, for setting domain limits in special functions, and for many other purposes in scientific computation, it is essential to know machine [Epsilon] (the floating-point precision) and machine [Omega] (the range of floating-point values). The landmark papers [Fox et al. 1978a; 1978b] gave a comprehensive and reliable solution that has since been widely adopted in the mathematical software community. They tabulated machine constants as defined by the Brown floating-point model [Brown 1977; 1981; Brown and Feldman 1980] for computers then in existence. The constants were expressed as Fortran data statements, and all were "commented out" by placing a "C" in column 1. Users were instructed to locate the particular set of constants for their local machine and to replace the "C" by a space before compilation.

These Fortran subroutines were a significant advance, not only because they quickly became a de facto standard, but because the machine epsilons tabulated there were based on extensive tests of each computer's floating-point functional units [Schryer 1984]. Not infrequently in that era, hardware (and, especially, double-precision software) was sloppy enough to require penalizing by several bits from what one would naively guess by counting bits in the floating-point representation.

The infelicity of this approach was the need for editing before compilation. As mathematical software and compilers continued to improve, these edits soon became the only changes needed for compilation and as such stood out as an annoyance. Also, the age of centralized computing with precompiled libraries administered by professional staff gave way to small computers with numerical software compiled as needed, and run by engineers with only occasional interest in mathematical software. Users voiced a strong preference for code that would compile correctly on heterogeneous machines without any change to source and without any compiler options.

So mathematical software started appearing again with the approach that had been used in earlier times: computing machine epsilon as needed by executing a loop like while (1 + [Epsilon] [is greater than] 1) {[Epsilon] = [Epsilon]/2}. This had several disadvantages. Efficiency was the least of these, since one could easily enough arrange to compute epsilon only once per program execution. Worse was the "arms race" between software writers and compiler optimizers, as attempts to store 1 + [Epsilon] in memory to avoid extra-precise registers led to global variables, leaf functions, and even more elaborate tricks to fool the increasingly clever compilers. Other problems were the inability to compute machine [Omega] and hidden assumptions in the loop calculation. (The most widely used "automatic" routine today, slamch in lapack [Anderson et al. 1992], says in its comments that it fails on Cyber machines because of assumptions about exponent limits being near a power of 2.) The worst disadvantage was not widely recognized: the simple loop fails to catch defects in the machine arithmetic that Schryer's more extensive floating-point tests find.

2. THE TRICK

We propose a solution that combines the automation general users desire and the authoritative results that experts require: look at bit patterns (in a portable way) to identify whether the machine is one of the known versions. Our implementation of r1mach and d1mach recognizes big- and little-endian IEEE, Vax D and G, IBM mainframe, Cray 1, T3E, and Convex C-1 formats, with possible compiler autodoubling. We have also updated i1mach, but observe that this function is less widely used. When our method does not identify the computer as one of the known types, it points to a separate file giving verified constants for machines believed to be no longer in use. The source for this new d1mach is 8KB, a noticeable savings over the previous implementation, which had grown to over 20KB for the slatec version or 35KB for the Algorithm 528 style with portable Fortan 66 output. This new implementation of the machine parameter functions has the same calling sequence as the previous implementation and should give bit-for-bit identical results.

The essence of the trick is this. Assign the value 1234567 to a real variable, and inspect an integer variable equivalenced to it. If the integer is 1234613304, then we have good evidence that the machine uses IEEE format. The other floating-point formats give rise to other bit patterns that can be recognized equally easily.

The main difficulty is to express the necessary 64-bit integers for Cray 1 and T3E in source code that will compile on machines with 32-bit integers. This could be accomplished with double-precision (floating-point) values, except that a few machines complain about double-precision constants being unsupported. So we use an auxiliary subroutine and small constants, attempting to avoid grief with compilers that zealously optimize 32-bit integer expressions.

Passing data between equivalenced integer and real variables is neither prohibited nor defined by the Fortran 77 standard, but all systems we have encountered allow it and implement the obvious sharing.

One could imagine a hypothetical machine that used IEEE floating-point representation, but to save transistors sometimes mangled the low-order bits of the mantissa. In the unlikely event that such a machine garners any market success, we will update our implementation to recognize it by testing for the offending operations. A more difficult case would be a system with heavily optimizing compilation on a front-end machine having different arithmetic than the main execution hardware; fortunately, such systems tend to be programmed in modern languages and by experts sensitive to porting issues.

We believe our approach provides an automatic and rigorously correct answer for all computers in wide use today, though "wide use" is perhaps a subjective term. One colleague has asserted that the Soviet BESM 6 is widely used but not supported by d1mach. For such machines, the former practice of editing data statements is retained.

3. HISTORICAL NOTE ON I1MACH (6)

In 1984, under the direction of Phyllis Fox and with contributions from various colleagues at Bell Labs, we completed the Third Edition of the PORT subroutine library, the core of which is Algorithm 528, the PORT frame. In doing so, we looked ahead to the coming availability of Fortran 77 compilers and adjusted the PORT frame so a simple preprocessor could change its source from Fortran 66 to Fortran 77 format. The main difference is that the Fortran 77 format uses character strings for error messages, rather than Hollerith data stored in integer variables. In making this change, we thought it reasonable to interpret i1mach(6) as the number of characters per character storage unit; for Fortran 66, this is the number of Hollerith characters that fit into an integer, and for Fortran 77 it is simply 1. Accordingly, for over 13 years we distributed source for an i1mach whose Fortran 77 version returned 1 for i1mach(6). This caused no trouble, probably because i1mach(6) is rarely used outside of a couple of routines in the PORT frame.

Unfortunately, other versions of i1mach, such as those distributed with the slatec library, have continued to interpret i1mach (6) as the number of chacters that fit into a Fortran integer. To avoid possible confusion, we recently adjusted our i1mach and the PORT frame routines (e9rint and seterr) affected by i1mach(6) so that the Fortran 77 versions of the frame routines no longer reference i1mach (6), and so i1mach (6) = 4 on IEEE arithmetic machines, where Fortran integers typically occupy four bytes.

4. THE FUTURE

The new versions of d1mach, r1mach, and i1mach have been deposited in netlib in the "blas" chapter, replacing the older versions there. The constants for older machines are in the file old1mach there.

It is worth emphasizing that the entire problem disappears in modern languages such as Ada or C or Fortran 9x. These languages provide built-in environmental enquiries, which should obviously be used. For quick conversion of old Fortran 77 programs to Fortran 9x, Bo Einarsson has created Fortran 90 versions using the language intrinsic functions.(1) Nevertheless, Fortran 77 is still in sufficiently widespread use that an improved d1mach is of more than historical interest.

ACKNOWLEDGMENTS

We thank Abhijit Bose, Kevin Cahill, Bo Einarsson, Fred Fritsch, Sven Hammarling, Norm Schryer, Layne Watson, and Warren Wiscombe for help testing versions of the code.

(1) Available via http://www.nsc.liu.se/~boein/ifip/kyoto/workshop-info /proceedings/einarsson/d1mach.html.

REFERENCES

ANDERSON, E., BAI, Z., BISCHOF, C., DEMMEL, J., DONGARRA, J., DO CROZ, J., GREENBAUM, A., HAMMARLING, S., MCKENNEY, A., OSTROUCHOV, S., AND SORENSEN, D. 1992. LAPACK User's Guide. SIAM, Philadelphia, PA.

BROWN, W. S. 1977. A realistic model of floating-point computation. In Mathematical Software III, J. R. Rice, Ed. Academic Press, Inc., New York, NY, 343-360.

BROWN, W. S. 1981. A simple but realistic model of floating-point computation. ACM Trans. Math. Softw. 7, 4 (Dec.), 445-480.

BROWN, W. S. AND FELDMAN, S. I. 1980. Environment parameters and basic functions for floating-point computation. ACM Trans. Math. Softw. 6, 4 (Dec.), 510-523.

Fox, P. A., HALL, A. D., AND SCHRYER, N. L. 1978a. Algorithm 528: Framework for a portable library. ACM Trans. Math. Softw. 4, 2 (June 1978), 177-188.

Fox, P. A., HALL, A. D., AND SCHRYER, N. L. 1978b. The PORT mathematical subroutine library. ACM Trans. Math. Softw. 4, 2 (June 1978), 104-126.

SCHRYER, N. L. 1984. Determination of correct floating-point model parameters. In Sources and Development of Mathematical Software, W. R. Cowell, Ed. Prentice-Hall Series in Computational Mathematics. Prentice-Hall, Englewood Cliffs, NJ, 360-366.

Received: August 1996; revised: July 1998 and December 1998; accepted: December 1998

This material was presented under the title "d1mach revisited: no more uncommenting data statements" at the IFIP WG2.5 International Workshop on "Current Directions in Numerical Software and High Performance Computing," October 19-20, 1995, Kyoto, Japan. Authors' address: Lucent Technologies, Bell Labs, Murray Hill, NJ 07974; email: ehg@research.bell-labs.com.

Permission to make digital/hard copy of part or all of this work for personal or classroom use is granted without fee provided that the copies are not made or distributed for profit or commercial advantage, the copyright notice, the title of the publication, and its date appear, and notice is given that copying is by permission of the ACM, Inc. To copy otherwise, to republish, to post on servers, or to redistribute to lists, requires prior specific permission and/or a fee.

Printer friendly Cite/link Email Feedback | |

Author: | GAY, DAVID M.; GROSSE, ERIC |
---|---|

Publication: | ACM Transactions on Mathematical Software |

Date: | Mar 1, 1999 |

Words: | 1762 |

Previous Article: | Remark on Algorithm 702--The Updated Truncated Newton Minimization Package. |

Next Article: | Complex Fans: A Representation for Vectors in Polar Form with Interval Attributes. |

Topics: |