Intel® Fortran Compiler
Build applications that can scale for the future with optimized code designed for Intel® Xeon® and compatible processors.
Announcements
Welcome to the Intel Community. If you get an answer you like, please mark it as an Accepted Solution to help others. Thank you!

## Floating point equivalence check

Beginner
232 Views

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....

13 Replies
Black Belt
232 Views

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```

Employee
232 Views

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.

Beginner
232 Views

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....

Beginner
232 Views

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  ```

Employee
232 Views

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".

Beginner
232 Views

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

Employee
232 Views

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.

Black Belt
232 Views

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

Employee
232 Views

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.

Black Belt
232 Views

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.

Black Belt
232 Views

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?

Black Belt
232 Views

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!

Black Belt
232 Views

>>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