Path: utzoo!attcan!uunet!ncrlnk!ncr-sd!hp-sdd!hplabs!ucbvax!CIM-VAX.HONEYWELL.COM!derstad
From: derstad@CIM-VAX.HONEYWELL.COM ("DAVE ERSTAD")
Newsgroups: comp.sys.apollo
Subject: Fortran 9.95 problem
Message-ID: <8811231951.AA29470@umix.cc.umich.edu>
Date: 23 Nov 88 18:50:00 GMT
Sender: daemon@ucbvax.BERKELEY.EDU
Organization: The Internet
Lines: 56


Here's one for the die-hard Fortran users out there.  Surely, someone
has run into this.

We're running a proprietary version of the UC Berkely SPICE program, 
a large Fortran program for analysis of integrated circuits.  The 
part we're having problems with is part of the public Domain code,
not our extensions, and has functioned properly on many previous 
releases of Aegis, the fortran compiler, etc., not to mention
virtually every other computer in the world.

With the 9.95 version of the compiler, we have problems.  Essentially,
there is a (actually many) segments of the code which look like

    A = something
    B = DSQRT(A*A + C*C)
    D = A + B - C
    IF D = 0 jump over some code

What is happening is that A will many times evaluate to identically
0.  B is then set to the same value as C without any problem.  However,
due to the way the optimized code works, the D evaluation is performed
with B in an 80 bit floating point register, but C comes from a 64
bit memory location.  (all variables are REAL*8).  Thus, D becomes 
non-zero.  Later, it gets used in an expression like
    E = 1.0D0 + D*something
    F = 1.0 / (1 - E)

which of course aborts with a divide by 0 error.

Now, I don't necessarily like the way this is coded, but it has worked
successfully at thousands of different sites on probably dozens of
different computer systems.  It devi
It definitely works if, with a given platform, DSQRT(C*C) is identically
C.  I think (but have not gone through and verified) that if it is
evaluated with a fixed precision throughout, whether it be 64 bits
or 80 bits, that it works as well.  Even if DSQRT(C*C) - C has 
round-off error, E becomes non-one and the code works.  What 
is causing the problem here is that DSQRT(C*C) - C is determined
to be non-zero based on 80 bit arithmetic, while the later
evaluation of E uses a 64 bit representation of DSQRT(C*C) - C
which rounds to zero.

I should mention that when I say the code "works" when E is non-one,
it may still yield screwy results.  However, the SPICE algorithms 
are capable of correcting for those screwy results.  Essentially,
SPICE makes some pretty strange guesses as to how the circuit is
operating, but is able to recover from them.

Anyone run into anything similar?  Any help would be appreciated.  So
far, the only definite workaround we know is to shut off cpu optimization;
however, then we lose much of the advantage of the hardware floating 
point support.

Dave Erstad
Honeywell SSED