Floating point equivalence check

[chauvjo]

I know this topic has been visit many times in this forum but I was wondering if there is now a standard way of coding floating point equivalence in an if statement. Consider the following:

if (diff == 0.0_dp) then
end if

Does the Fortran standard (Fortran 2008 and below) now include a formal way of coding this evaluation?  Maybe as special version of the if statement:

ifsame (diff,0.0_dp) then
end if

If not, what is "best practice" approach to this problem?

Thanks....

[IanH]

There is no "standard" way.

Given Fortran real numbers are an approximation of really real numbers, what do you (given your problem domain and algorithm) consider to be the same as zero (or "no difference")? 

Figure that out, then write the comparison appropriately.  It may help to put the test in a function.

IF (close_enough_to_zero(diff)) THEN
  ...


PURE FUNCTION close_enough_to_zero(x)
  REAL(dp), INTENT(IN) :: x
  LOGICAL :: close_enough_to_zero
  
  ! I'm not particularly fussy.
  close_enough_to_zero = ABS(x) < 1E6_dp

[Steven_L_Intel1]

Keep in mind that if the values are large, even a one-bit difference in the fraction might be bigger than some fixed value. (I hope Ian meant 1E-6 instead of 1E6!) Standard Fortran has intrinsics such as SPACING and RRSPACING that can help you find tolerance values that are appropriate for the magnitude of the values. You might want to pass the actual values to a function rather than the difference if this is a concern for you.

[chauvjo]

Would it be possible to see an example that would handle floating point comparison of any two floating point values for equivalence not just to zero?    Thanks....

[chauvjo]

I just found the following in one of our codes.  Any thoughts or suggestions or modifications?

      LOGICAL FUNCTION dpc(X,Y)      ! Double Precision Comparison Function

C--------------------------------------------------------------------------------
C   ULP:         Unit of data precision. The acronym stands for "unit in 
C                the last place," the smallest possible increment or decrement 
C                that can be made using a machine's floating point arithmetic. 
C                A 0.5 ulp maximum error is the best you could hope for, since 
C                this corresponds to always rounding to the nearest representable 
C                floating-point number. Value must be positive - if a negative 
C                negative value is supplied, the absolute value is used. 
C                If not specified, the default value is 1. 
C--------------------------------------------------------------------------------
      
      IMPLICIT NONE
C ..
C .. Scalar Arguments ..
      DOUBLE PRECISION X, Y
C ..
C .. Parameters ..
      DOUBLE PRECISION ulps
      PARAMETER(ulps=5.0D0)     ! Tolerance on comparison of two real variables 
C..
C .. Executable Statements ..

      dpc = ABS(X-Y) < ulps*SPACING(MAX(ABS(X),ABS(Y)))         ! eps = ulps*SPACING(MAX(ABS(X),ABS(Y)))
      
      RETURN
      END FUNCTION dpc  

[Steven_L_Intel1]

Yeah - that's the way to do it - get your ulps for nothing and your diffs for free. Indeed, that's a very fine method. It uses SPACING to determine the largest spacing between two representable values (ULPs) and multiplies by the ULPs it considers "close enough".

[chauvjo]

Thanks Steve...How do I know if the ulps value (which is 5 in the sample code) is correct?

[Steven_L_Intel1]

There is no "correct" here - it's up to you to decide what is "close enough", maybe experiment with some of your real data. A lot depends on the precision and accuracy of your input values and what sort of transformations are done to them. 5ULPs strikes me as an eminently reasonable starting point.

[jimdempseyatthecove]

Why doesn't the documentation for SPACING and EPSILON have "See Also" hyperlinks to each other.

SPACING(X) ... ABS(X)*EPSILON(X)

Though SPACING should be more efficient as it requires only an AND and INC (or OR with 1)

Jim Dempsey

[Steven_L_Intel1]

I'll make the suggestion that EPSILON, SPACING and RRSPACING link to each other. But we don't want to get carried away with this.

[IanH]

Steve Lionel (Intel) wrote:
(I hope Ian meant 1E-6 instead of 1E6!)

As I said, I'm not fussy.  Plus my models all converge real quick.

[mecej4]

ianh wrote:

Quote:

Steve Lionel (Intel) wrote:

(I hope Ian meant 1E-6 instead of 1E6!)

As I said, I'm not fussy.  Plus my models all converge real quick.

With such a relaxed criterion, they (your models) are already "there", so who needs algorithms?

[IanH]

mecej4 wrote:

Quote:

ianh wrote:

Quote:

Steve Lionel (Intel) wrote:
(I hope Ian meant 1E-6 instead of 1E6!)

As I said, I'm not fussy.  Plus my models all converge real quick.

With such a relaxed criterion, they (your models) are already "there", so who needs algorithms?

I put that sort of superfluous complexity aside long ago. 

My approach means that initial guesses are pretty important though - but that has its benefits as a consultant... I just need to get the client to remind me of the answer that they wanted the model to give.

Time to get back to writing my new accounting software!

[jimdempseyatthecove]

>>I just need to get the client to remind me of the answer that they wanted the model to give.

Time to get back to writing my new accounting software!

Don't you mean the two sets of answers they want ?)

IOW Accounting with complex numbers: one real, the other imaginary

Jim Dempsey