(1.0+eps) .ne. 1.0

Nelson H. F. Beebe sun!science.utah.edu!Beebe
Mon Jun 11 09:50:20 PDT 1990


The use of an expression like

(1.0+eps) .ne. 1.0

to test for machine epsilon, or convergence of an iteration,
raised the question "What if the compiler evaluates the
expression in higher precision, such as IEEE temporary
real?".

This indeed happens on more than one compiler I have used on
IEEE arithmetic systems (notably Intel 808x7 processors); it
has also happened on older architectures, such as Honeywell,
that had registers longer than the word size.

The LINPACK User's Guide, written in 1978, says (I.4):

>> The convergence of the iteration in the singular value
>> decomposition is tested in a machine-independent manner
>> by statements of the form
>>     TEST1 = something not small
>>     TEST2 = TEST1 + something possibly small
>>     IF (TEST1 .NE. TEST2) ...

The actual code from SSVDC has

>>          TEST = ABS(S(L)) + ABS(S(L+1))
>>          ZTEST = TEST + ABS(E(L))
>>          IF (ZTEST .NE. TEST) GO TO 380
>> ...
>>          ZTEST = TEST + ABS(S(LS))
>>          IF (ZTEST .NE. TEST) GO TO 420

I couldn't find other such instances in the rest of
LINPACK by searching for .EQ. and .NE..  

EISPACK/2 has over a dozen routines in which explicit
assignments of a machine-dependent MACHEP or RADIX occur; I
changed our local sources a decade ago to call subroutines
to get these values.

The workaround for systems that compute intermediate
expressions in higher precision (which is almost always a
good thing anyway) is to write

     TEST1 = something not small
     TEST2 = TEST1 + something possibly small
     CALL STORE (TEST1)
     CALL STORE (TEST2)
     IF (TEST1 .NE. TEST2) ...

where STORE is a SEPARATELY-COMPILED routine that can in
fact do nothing to its argument; the compiler is forced to
assume that TEST1 and TEST2 may have been changed by the
calls, and therefore, to have stored them in memory, instead
of preserving them in longer registers.  The IF (...) then
performs as it would on a machine in which registers and
memory words have the same sizes.

Some compilers, such as Lahey F77L on the IBM PC, offer a
compile-time option to force memory storage of intermediate
results; while that is helpful in cases like the above, it
is also prone to being forgotten, especially in code that is
distributed to other users.  Thus, a source-level change
such as CALL STORE() seems to be the best defensive
approach.
-------



More information about the Numeric-interest mailing list