Given the nice inquiry functions (EPSILON, TINY, etc.) provided by modern Fortran, one may not appreciate how much effort was needed to adapt a Fortran 66 or Fortran 77 program to a new computer, or to port a program from one architecture to another. Today, I came across a paper from Y2K-1 that addressed the question, “Does my computer use IEEE?” The F77 program below implements the trick given in the paper to answer the question. The last line contains a link to the paper.
program AmIEEE
implicit none
integer i
real r
equivalence (i,r)
r = 1234567
print 10, i .eq. 1234613304
10 format(' Does this machine use IEEE FP? ',L7)
end
!Self-Adapting Fortran 77 Machine Constants: Comment on Algorithm 528
! by David M. Gay and Eric Grosse
!ACM TOMS Vol. 25, No. 1, March 1999, Pages 123–126
!https://dl.acm.org/doi/pdf/10.1145/305658.305711
P.S.: Edited in response to comments below. “Implicit None” was recognized by most F77 compilers after persuasion by the military ( MIL-STD 1753 ) and lowercase letters were tolerated (see an interesting description of the extensions and features, in which lowercase letters are named “almost standard”).
Thanks for identifying this paper. It is very interesting.
I never really had much use for this detailed approach, as I have limited my calculations to real8 and the occasional accumulator as real10 (which I assumed complied with IEEE 754). Typically I just hoped for the best precision available, while itterative convergence tests were typically for much less precision based on the computation algorithm, than the precision of the numeric format.
Do you have any links to examples of “slamch in lapack” or “updated i1mach,” or “this new d1mach” as referred to in the paper at 1998 or the new versions of “d1mach, r1mach, and i1mach h” ?
It would be interesting to understand the history of these constants, although my impression was these became less significant after IEEE Standard 754 and the release of Fortran 90.
How does all this relate to Ifort not fully supporting ISO/IEC/IEEE 60559:2011, if I understand this correctly ?
@JohnCampbell ,
The “almost portable” versions of i1mach(), r1mach() and d1mach() may be found in the Netlib directory port . The codes there have special-cased Honeywell DPS, PDP 11 and Univac 1100 machines, and all the rest are treated using portable code. The earlier versions had the feature of the user having to remove ‘C’ in the first column of the lines that were relevant to the user’s machine. We still run into those when we download old Fortran code from, for example, TOMS. If interested in the historical perspective, download, adapt and run TOMS-528 and TOMS-665.
The Lapack auxiliary routines slamch() and dlamch() use a single character argument instead of an integer, but these routines appear to be rarely used in published source codes, in my experience.
I work with old Fortran software quite often. I have found it convenient to throw out any versions of ?1mach() routines in the downloaded files, and use the carefree versions below.
double precision function d1mach (i)
integer i
if (i .lt. 1 .or. i .gt. 5) call xermsg ('SLATEC', 'D1MACH',
+ 'I OUT OF BOUNDS', 1, 2)
select case (i)
case (1)
d1mach = tiny(1d0)
case (2)
d1mach = huge(1d0)
case (3)
d1mach = epsilon(1d0)/2
case (4)
d1mach = epsilon(1d0)
case (5)
d1mach = log10(2d0)
end select
return
end
A copy of the above, with “real”, “r1” and “e0” replacing “double precision”, “d1” and “d0”, will constitute r1mach().
Here is i1mach():
integer function i1mach (i)
integer i, imach(16), output
save imach
equivalence (imach(2),output)
data imach( 1) / 5 /
data imach( 2) / 6 /
data imach( 3) / 0 /
data imach( 4) / 0 /
data imach( 5) / 32 /
data imach( 6) / 4 /
data imach( 7) / 2 /
data imach( 8) / 31 /
data imach( 9) / 2147483647 /
data imach(10) / 2 /
data imach(11) / 24 /
data imach(12) / -125 /
data imach(13) / 128/
data imach(14) / 53 /
data imach(15) / -1021 /
data imach(16) / 1024/
if (i .lt. 1 .or. i .gt. 16) go to 10
i1mach = imach(i)
return
10 continue
write (unit = output, fmt = 9000)
9000 format ('*** ERROR 1 IN I1MACH - I OUT OF BOUNDS')
stop
end
I notice that this article does not use the tiny(), epsilon(), or spacing() intrinsic functions. nearest() is another very useful function for writing portable code, and it has no r1mach() or d1mach() equivalent. All of these fortran intrinsic functions work for any real kind, not just the old single and double precision. Also, just in case it might make a difference, the intrinsics can be used in parameter statements at compile time (and they are likely evaluated at compile time even within expressions), whereas r1mach() and d1mach() are evaluated at run time.
This is a much better way to specify these constants with modern fortran. As for programming style, no one should be using d exponents these days to specify the kind value, just use the flexible kind facility of fortran. This would allow these functions to be written within a generic interface that allows all real kinds, not just the old single and double precision. Also, it might be a good idea to declare the functions pure, as this might facilitate compiler optimizations in expressions that reference these functions.
As I said in a separate post, there are also advantages to just using the fortran intrinsics directly in codes, and avoid the r1mach/d1mach run time overhead entirely.
Note the authors of the paper cited in the original post wrote back in 1999:
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.
I suppose OP might be trying to convey the statement, “Nevertheless, Fortran 77 (sic) is still in sufficiently widespread use that an improved d1mach is of more than historical interest,” remains relevant in year 2023?
Yes I agree. I typically don’t use the *1mach routine in new code, or when I modernize old code. I just replace the calls with parameters. Then if I change the real kind it all “just works”.