Skip to content

Add a precision argument to mo_utils comparison operators

I wrote a little program, simulating the behaviour of == and eq from mu_utils.f90 for a set of real numbers -huge(dp), -1.0, 0.0, 1.0, huge(dp) altered by -1.0e-9_dp, -epsilon(dp), epsilon(dp), 1.0e-9_dp.

module mo_kind
  implicit none
contains
  function eq(a, b) result(boolean)
   use iso_fortran_env, only: sp=>real32, dp=>real64, qp=>real128
   use iso_fortran_env, only: i1=>int8, i2=>int16, i4=>int32, i8=>int64
   
   real(dp), intent(in) :: a, b
   logical :: boolean
   
   if ((epsilon(1.0_dp) * abs(b) - abs(a - b)) < 0.0_dp) then
      boolean = .false.
   else
      boolean = .true.
   end if
 end function eq
end module mo_kind

 program float_comparison
   use iso_fortran_env, only: sp=>real32, dp=>real64, qp=>real128
   use iso_fortran_env, only: i1=>int8, i2=>int16, i4=>int32, i8=>int64
   use mo_kind, only: eq

   implicit none
   real(dp), parameter :: as(5) = [huge(1.0_dp) * (-1.0_dp) + 10.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, huge(1.0_dp) - 10.0_dp]
   real(dp), parameter :: mods(4) = [-1.0e-9_dp, -epsilon(1.0_dp), epsilon(1.0_dp), 1.0e-9_dp]
   real(dp) :: a, b
   integer(i4) :: i, j
   logical :: boolean_1, boolean_2

   do i=1, size(as)
      a = as(i)
      do j=1, size(mods)
         b = a + mods(j)
         print*, 'mo_utils: eq(', a, ', ', b, ') = ', eq(a,b)
         print*, 'mo_utils: eq(', b, ', ', a, ') = ', eq(b,a)
         print*, '==      :   ', a, '== ', b, ' = ', a == b
      end do
   end do
 end program float_comparison

which generates this output (GFortran9.3 on MacOSX)

 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 ==      :     -1.7976931348623157E+308 ==   -1.7976931348623157E+308  =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 ==      :     -1.7976931348623157E+308 ==   -1.7976931348623157E+308  =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 ==      :     -1.7976931348623157E+308 ==   -1.7976931348623157E+308  =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 mo_utils: eq(  -1.7976931348623157E+308 ,   -1.7976931348623157E+308 ) =  T
 ==      :     -1.7976931348623157E+308 ==   -1.7976931348623157E+308  =  T
 mo_utils: eq(  -1.0000000000000000      ,   -1.0000000010000001      ) =  F
 mo_utils: eq(  -1.0000000010000001      ,   -1.0000000000000000      ) =  F
 ==      :     -1.0000000000000000      ==   -1.0000000010000001       =  F
 mo_utils: eq(  -1.0000000000000000      ,   -1.0000000000000002      ) =  T
 mo_utils: eq(  -1.0000000000000002      ,   -1.0000000000000000      ) =  T
 ==      :     -1.0000000000000000      ==   -1.0000000000000002       =  F
 mo_utils: eq(  -1.0000000000000000      ,  -0.99999999999999978      ) =  F
 mo_utils: eq( -0.99999999999999978      ,   -1.0000000000000000      ) =  T
 ==      :     -1.0000000000000000      ==  -0.99999999999999978       =  F
 mo_utils: eq(  -1.0000000000000000      ,  -0.99999999900000003      ) =  F
 mo_utils: eq( -0.99999999900000003      ,   -1.0000000000000000      ) =  F
 ==      :     -1.0000000000000000      ==  -0.99999999900000003       =  F
 mo_utils: eq(   0.0000000000000000      ,   -1.0000000000000001E-009 ) =  F
 mo_utils: eq(  -1.0000000000000001E-009 ,    0.0000000000000000      ) =  F
 ==      :      0.0000000000000000      ==   -1.0000000000000001E-009  =  F
 mo_utils: eq(   0.0000000000000000      ,   -2.2204460492503131E-016 ) =  F
 mo_utils: eq(  -2.2204460492503131E-016 ,    0.0000000000000000      ) =  F
 ==      :      0.0000000000000000      ==   -2.2204460492503131E-016  =  F
 mo_utils: eq(   0.0000000000000000      ,    2.2204460492503131E-016 ) =  F
 mo_utils: eq(   2.2204460492503131E-016 ,    0.0000000000000000      ) =  F
 ==      :      0.0000000000000000      ==    2.2204460492503131E-016  =  F
 mo_utils: eq(   0.0000000000000000      ,    1.0000000000000001E-009 ) =  F
 mo_utils: eq(   1.0000000000000001E-009 ,    0.0000000000000000      ) =  F
 ==      :      0.0000000000000000      ==    1.0000000000000001E-009  =  F
 mo_utils: eq(   1.0000000000000000      ,   0.99999999900000003      ) =  F
 mo_utils: eq(  0.99999999900000003      ,    1.0000000000000000      ) =  F
 ==      :      1.0000000000000000      ==   0.99999999900000003       =  F
 mo_utils: eq(   1.0000000000000000      ,   0.99999999999999978      ) =  F
 mo_utils: eq(  0.99999999999999978      ,    1.0000000000000000      ) =  T
 ==      :      1.0000000000000000      ==   0.99999999999999978       =  F
 mo_utils: eq(   1.0000000000000000      ,    1.0000000000000002      ) =  T
 mo_utils: eq(   1.0000000000000002      ,    1.0000000000000000      ) =  T
 ==      :      1.0000000000000000      ==    1.0000000000000002       =  F
 mo_utils: eq(   1.0000000000000000      ,    1.0000000010000001      ) =  F
 mo_utils: eq(   1.0000000010000001      ,    1.0000000000000000      ) =  F
 ==      :      1.0000000000000000      ==    1.0000000010000001       =  F
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 ==      :      1.7976931348623157E+308 ==    1.7976931348623157E+308  =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 ==      :      1.7976931348623157E+308 ==    1.7976931348623157E+308  =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 ==      :      1.7976931348623157E+308 ==    1.7976931348623157E+308  =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 mo_utils: eq(   1.7976931348623157E+308 ,    1.7976931348623157E+308 ) =  T
 ==      :      1.7976931348623157E+308 ==    1.7976931348623157E+308  =  T

I learned some things:

  1. a == a is always True, trivial but important
  2. eq(a, a + epsilon) != eq(a + epsilon, a) at least for the (edge) cases 1.0_dp and -1.0_dp
  3. a == a + 1e-9 is True when a tends towards huge(dp)

I think this leads to the following:

  1. the usage of eq from mo_utils.f90 is not safe in all circumstances
  2. for the common case, where you init some value (e.g. a = nodata_dp) and later check for it with if (a == nodata_dp), it is okay to use == instead of eq. This of course does not hold true, if you read values a from a file because of precision loss.
  3. We need to improve the eq functions in mo_utils.f90 to check for absolute or relative tolerance (defaulting to epsilon in order to keep backwards compatibility. So something like abs(a - b) <= (atol + rtol * absolute(b)) (taken from numpy.allclose).

What do you think @muellese , @thober, @shresthp , @kaluza?

Edited by Sebastian Müller