forums.silverfrost.com

eric_carwardine

G'day, folks Smile

What's a more reliable way of comparing floating point quantities for inequality? The ftn95 compiler reports that -

davidb · Joined: 17 Jul 2009 Posts: 560 Location: UK

LitusSaxonicum · Posted: Sun May 08, 2011 7:13 pm Post subject:

Hi Eric,

This warning often irritates me too, particularly as I try to catch divide-by-zeroes that are the result of wrong data input (normally it hasn't been done) and not a calculation I expect to produce a "perfect" answer - I've known about roundoff for more than 4 decades, and already know that in some cases these tests aren't a good idea!

If you set a value in an assignment statement, or read a value in from a datafile, and the value is zero (or any other value for that matter), then it is perfectly proper to test for it.

e.g.

JohnCampbell · Joined: 16 Feb 2006 Posts: 2615 Location: Sydney

Eric,

I have tried for a while to adopt an approach which limits the warnings that FTN95 provide.
I have 2 cases where I adopt different approaches.

1) The first, where I choose to ignore the error, such as when:

A = B

{do some calculations}

if (A == B) ... ! check if A has not changed

In this case B is a special value which I have used to previously apply to A ( B can also be 0 ). I am testing if A has changed in any way. In this case I ignore the warning.
I've also asumed that "zero" is a unique real number.

2) Where the warning has some meaning, where you are testing the value of A and must recognise that real arithmetic has some finite accuracy. In this case you need to find an error "ACC" which is within limits of the calculation or within the accuracy of the method. Examples are:
1.0d-12 or 1.0e-6 may be realistic for real*8 or real*4, but also depend on the magnitude of A, as if A is expected to be 1.e-6, then a test for 1.e-16 may be more sensible for real*8 calculations.
In short it is difficult to apply a general rule for ACC.

For Eddie's example
IF (DABS(A-B) .GT. ACC) THEN
A sensible approach to ACC could be ACC = ( Abs(A) + Abs(b) ) *1.d-8,

It all depends on how frequently you do the calculation and what is the consequence of getting the test wrong. In the case of avoiding divide by zero for IF (A.NE.B) is true and if A-B is very small, and 1 / (A-B) will be very large. You need to check the consequence, which is often not a problem.

Other problems come in when ACC is set to a rounded value,
eg ACC = 0.1 .
In this case 0.1_1 is different to 0.1_2 and this can cause some other problems. Try the following code:

Wilfried Linder · Posted: Mon May 09, 2011 9:08 am Post subject:

Only as an additional remark: In some cases the tiny function may be helpful which gives the smallest positive number bigger than zero. This is really near to zero, so I usually take something like r = 10.D0*tiny(1.D0) and then the expression "if (abs(x) .lt. r) ..."

Regards - Wilfried

eric_carwardine · Posted: Wed Jun 01, 2011 2:21 pm Post subject:

G'day, folks Smile

Crikey, it's June already! I had intended getting back here much sooner. Many thanks for all the comments!

davidb wrote:
" ... and hope the compiler doesn't spot what you are doing ... "
Almost sounds like we're dealing with HAL, David, from the film "2001: A Space Odyssey" Smile

I wonder if the FTN95 compiler will ever pass the Turing Test?

Eddie wrote:
" ... While programming with Clearwin I have just about given up on real*4 ... "
I must admit, Eddie, that Clearwin's boycotting of real*4 caused me to frown a bit ... well, a lot. Then I wondered if I could get around it by using Fortran's 'internal write'. The essential code would be ...