Re: CDC Floating point = tests?

stocksami@xxxxxxxxxxxxx wrote:
On Jan 30, 8:54 pm, nos...@xxxxxxxxxxxxx (Richard Maine) wrote:
stocks...@xxxxxxxxxxxxx <stocks...@xxxxxxxxxxxxx> wrote:
On the old CDC machines running EKS you could use a regular if test to
see if two floating point values were equal, and it worked! I've
ported lots of code from EKS to other platforms and I've had to change
those lines of codes many many times.
I'm not sure exactly what your question really is. First, it sounds as
though you are surprised that testing floating point values for equality
can ever work. Well, it can. No special compiler tricks are necessary.
It is a simple operation that any Fortran compiler can do - and can do
correctly. I half wonder whether you might be seeing warning messages
from other compilers and overreacting to them. There are compilers that
will give you a warning for every attempt to compare two floating point
numbers for equality. That is only a warning, suggesting that this might
be a bad idea. It does *NOT* mean that it won't work; it will work.
Whether it is a good idea is a separate matter. Sometimes it is a
perfectly fine idea. Other times it is not.

The other posters have concentrated on the question of why it is
sometimes not a good idea, what can be done about those cases, why the
results of the problematic cases can vary among compilers, and related
matters.

But your question seemed to me to have a hint of suprise that any code
could ever work with floating point equality tests. Perhaps I misread
that hint, but I thought it worth clarifying the sim ple case first. I
have seen people who had perfectly valid and "safe" code, but thought
that they had to change it because the compiler gave them a warning
message. The people did not understand the distinction between a warning
and an error. That might not be the case for you, in which case I
apologize for making a "trivial" point. But it seemed worth making just
in case.


My original question was prompted by results I have gotten on other
platforms, mostly AIX. I have never gotten any warnings, but the very
same code and input files that worked on EKS would fail on AIX because
the IF tests weren't reliable. I investigated and found that on AIX a
number such as 101.1 might have all sorts of internal representations
such as 101.099999999999999999 or 101.100000000000000001 resulting in
erratic results. Previous posts have mentioned that such tests can be
reliable if all of the library functions work similarly in how they
store numbers. The numbers were read in using formatted reads. One
set was from a file and another from internal reads of character
strings. There were also some floating point values that were
hardcoded, such as 0.0.

To get around the problem I replaced all such if tests with a function
call that returned true or false. The function would use a tolerance
value to determine equality.

For those cases where figuring a tolerance value is a hassle, have a lookee at what I use:

ELEMENTAL FUNCTION Compare_Real_Single( x, y, ulp ) RESULT( Compare )
REAL(Single), INTENT(IN) :: x
REAL(Single), INTENT(IN) :: y
INTEGER, OPTIONAL, INTENT(IN) :: ulp
LOGICAL :: Compare
REAL(Single) :: Rel
Rel = 1.0_Single
IF ( PRESENT( ulp ) ) THEN
Rel = REAL( ABS(ulp), Single )
END IF
Compare = ABS(x-y) < ( Rel * SPACING( MAX(ABS(x),ABS(y)) ) )
END FUNCTION Compare_Real_Single

with similar for double precision real and complex equivalents. All overloaded in a module like so:

PRIVATE
PUBLIC :: Compare_Float
INTERFACE Compare_Float
MODULE PROCEDURE Compare_Real_Single
MODULE PROCEDURE Compare_Real_Double
MODULE PROCEDURE Compare_Complex_Single
MODULE PROCEDURE Compare_Complex_Double
END INTERFACE Compare_Float

so that one just needs to do:

IF ( Compare_Float(x,y,ulp=ulp) ) THEN
....

for scalars, and something like

IF ( ALL(Compare_Float(x,y,ulp=ulp)) ) THEN
....

for arrays.

Works well.

BTW, from my header docs:

! James Van Buskirk and James Giles suggested this method for floating
! point comparisons in the comp.lang.fortran newsgroup.

just in case you think I came up with it. :o)

cheers,

paulv

--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
Ph: (301)763-8000 x7748
Fax:(301)763-8545
.

Relevant Pages

  • Re: CDC Floating point = tests?
    ... though you are surprised that testing floating point values for equality ... No special compiler tricks are necessary. ... I half wonder whether you might be seeing warning messages ...
    (comp.lang.fortran)
  • Re: CDC Floating point = tests?
    ... though you are surprised that testing floating point values for equality ... No special compiler tricks are necessary. ... I half wonder whether you might be seeing warning messages ...
    (comp.lang.fortran)
  • Re: Im not the only one who values interval-arithmetic
    ... I see rounding mode discussed in regard to floating ... It's only used by the compiler ... name of the package and something about how it works? ... coerced to an interval by mapping NUM -> before some ...
    (comp.lang.lisp)
  • Re: Hercules
    ... Although D-format floating point does use 64 bits to represent a floating point number, I don't call the instructions which process these numbers "64-bit instructions". ... In PL/I terms, instruction I received way back when indicated that FLOAT DECwas the way to go to use E-format in the compiled code, and FLOAT DECwas the way to go to use D-format in the compiled code. ... My reference to Fujitsu was simply to record my experience that Fj had gone to the trouble of making their PL/I compiler sufficiently "plug compatible" with IBM's Optimizing compiler that ...
    (comp.lang.pl1)
  • Re: Simple Random Number Generator/Function - Cliff RNG
    ... Tim and Alex: ... This has certainly informed me greatly on this 'PRNG'. ... PRNGs with the extra baggage of this floating point fuzziness! ... not just the compiler. ...
    (sci.math)