When 1.000000 is greater than 1.000000

[michael_green]

Hi All,

In the following code:

	            if(h2-h1 > defects(DefectIndex)%MaxLen)then
		            call DisplayError(hDlg,1,t,b,d,1,ErrorConfirmed)
	            end if

h2-h1 has the value 1.000000 and so does defects(DefectIndex)%MaxLen. Why then does the total expression evaluate to .TRUE. and why is my error message triggered?

h2 and h1 are locally defined Real(4), and defects(DefectIndex)%MaxLen (part of a global type) is also Real(4) whose actual value is read from a text file at run time.

This will be something about the way real numbers are encoded I expect, but how do I get around it?

Many thanks in advance,

Mike

[jimdempseyatthecove]

The numbers are not equal. There is a tiny binary difference.

As an example, x = 0.1, is not exactly representable in binary floating point. It would require infinite precision, and then would not be exact.

Do not be confused with what you see in the Debug watch window or with write statement because both of these will incorporate rounding. Display in hex and/or use Z format control to find what is actually use in your evaluations.

The decimal display of x=0.1 is 0.1000000, the Hex is #3DCCCCCD, where the exact value would replace the trailing "D" with an infinite number of "C"'s.

As for the fix, we know not the requirements of your code. As written, your code may or may not be giving you a false positive.

Jim Dempsey

[mecej4]

You can only "get around it" by modifying your expectations which, by the way, you have not told us. Remember that REAL variables are stored in a biased exponent, base-2 floating point representation. Real numbers such as 0.1_ten may have an infinitely long representation in base-2, and thus truncation errors can arise. You may find the following test program instructive.

program round1
implicit none
real :: x, h1=9871378.9, h2=9871377.8
integer :: i
character(len=3) :: str = '1.0'
!
read(str,*)x
if(x.eq.h1-h2)then
   print *, 'equal'
else
   print *, 'unequal'
endif
h1=1234.56
h2=12345.6
i=h2/h1
write(*,*)i
end

[Greg_T_]

Hi Mike,

Would you be able to add a small epsilon value within the greater-than comparison?  I've had success either subtracting a small epsilon to the quantity on the left, or adding a small epsilon to the quantity on the right of the greater-than comparison to avoid the inexact floating point value described above.  If you expect the values to have a large range of magnitudes, from very small to very large values, the epsilon value could be updated before the comparison by making it a small percentage tolerance of the values being compared.

Regards,
Greg

[jimdempseyatthecove]

As an alternate, if you are able to substitute different (smaller) units of measure, such that the majority of your numbers remain integral, then you may avoid many round-off errors (approximations). As an example, in the financial industry change your currency of dollars into cents or mills. In Physics, change meter to nanometer. The choice of unit might be critical, so you may need some experimentation. 

A similar situation occurs in summations of large numbers of values of large range of magnitudes. To mitigate loss of precision you would perform summations (reductions) on similarly ranged numbers into partial sums, then sum the partial sums. Alternately, you can sort the data from small to large then sum from small to large magnitude.

You can also change the way you calculate to produce more accurate results. Assume you want to find all values less than 10% of X

IF(A .LT. X / 10.0) ...
IF(A*10.0 .LT. X) ...

(assuming the compiler optimizations did not alter your expressions)
The second statement can be more precise than the first depending on the values.

Jim Dempsey

[van_der_merwe__ben]

We have the same sort of issue in our code and we handle it by adding small values when comparing two real numbers in some cases. There is generically no other way to handle the fact that numbers in practice are sometimes not exact.

We once had code which effectively came down to this:

flow1 = 0.9 * flowl + 0.1*flowv
flow2 = 0.1 * flowl + 0.9*flowv
if flow1 + flow2 < flowl + flowv then ERROR

which is just not good robust code, and would hit on some machines and not others, for some cases, but not all. Better to write:

if flow1 + flow2 + delta < flowl + flowv then ERROR

[Steven_L_Intel1]

That doesn't always work - see the example in the presentation I attached to your other thread. Better is to compare ABS of the difference to some epsilon value.