floating point equality wip

autotest/TestGenericUtils.f90

@@ -1,7 +1,7 @@
 module TestGenericUtils
   use testdrive, only: error_type, unittest_type, new_unittest, check
   use KindModule, only: DP
-  use GenericUtilitiesModule, only: is_same
+  use GenericUtilitiesModule, only: is_close
 
   implicit none
   private
@@ -12,41 +12,41 @@ module TestGenericUtils
   subroutine collect_genericutils(testsuite)
     type(unittest_type), allocatable, intent(out) :: testsuite(:)
     testsuite = [ &
-                new_unittest("is_same", test_is_same), &
-                new_unittest("is_same_both_near_0", test_is_same_both_near_0, &
+                new_unittest("is_close", test_is_close), &
+                new_unittest("is_close_both_near_0", test_is_close_both_near_0, &
                              should_fail=.true.), & ! expect failure for now, see below
-                new_unittest("is_not_same", test_is_not_same) &
+                new_unittest("is_not_close", test_is_not_same) &
                 ]
   end subroutine collect_genericutils
 
-  subroutine test_is_same(error)
+  subroutine test_is_close(error)
     type(error_type), allocatable, intent(out) :: error
 
     ! exact
-    call check(error, is_same(0.0_DP, 0.0_DP))
+    call check(error, is_close(0.0_DP, 0.0_DP))
     if (allocated(error)) return
 
     ! inexact (within tolerance)
-    call check(error, is_same(1.0000_DP, 1.0001_DP, eps=0.01_DP))
+    call check(error, is_close(1.0000_DP, 1.0001_DP, eps=0.01_DP))
     if (allocated(error)) return
-  end subroutine test_is_same
+  end subroutine test_is_close
 
-  subroutine test_is_same_both_near_0(error)
+  subroutine test_is_close_both_near_0(error)
     type(error_type), allocatable, intent(out) :: error
 
     ! relative comparison mode fails when a and b are close to 0
-    call check(error, is_same(0.0000_DP, 0.0001_DP, eps=0.01_DP))
+    call check(error, is_close(0.0000_DP, 0.0001_DP, eps=0.01_DP))
     if (allocated(error)) return
-  end subroutine test_is_same_both_near_0
+  end subroutine test_is_close_both_near_0
 
   subroutine test_is_not_same(error)
     type(error_type), allocatable, intent(out) :: error
 
-    call check(error, (.not. (is_same(1.0_DP, 1.0001_DP))))
+    call check(error, (.not. (is_close(1.0_DP, 1.0001_DP))))
     if (allocated(error)) return
 
     ! with tolerance
-    call check(error, (.not. is_same(1.0_DP, 1.0001_DP, eps=0.00005_DP)))
+    call check(error, (.not. is_close(1.0_DP, 1.0001_DP, eps=0.00005_DP)))
     if (allocated(error)) return
   end subroutine test_is_not_same
 

src/Model/Connection/GridSorting.f90

@@ -2,7 +2,7 @@ module GridSorting
   use KindModule, only: I4B, DP, LGP
   use ConstantsModule, only: DHALF
   use CellWithNbrsModule, only: GlobalCellType
-  use GenericUtilitiesModule, only: is_same
+  use GenericUtilitiesModule, only: is_close
   use BaseDisModule, only: dis_transform_xy
   implicit none
   private
@@ -75,11 +75,11 @@ function lessThan(n, m) result(isLess)
                     dis_bot_m(gcm%index))
 
       ! compare
-      if (.not. is_same(zn, zm, 10 * epsilon(zn))) then
+      if (.not. is_close(zn, zm, 10 * epsilon(zn))) then
         isLess = zn > zm
-      else if (.not. is_same(yn, ym, 10 * epsilon(yn))) then
+      else if (.not. is_close(yn, ym, 10 * epsilon(yn))) then
         isLess = yn > ym
-      else if (.not. is_same(xn, xm, 10 * epsilon(xn))) then
+      else if (.not. is_close(xn, xm, 10 * epsilon(xn))) then
         isLess = xn < xm
       else
         isLess = .false.

src/Model/GroundWaterFlow/gwf3csub8.f90

@@ -18,7 +18,7 @@ module GwfCsubModule
                              TABLEFT, TABCENTER, TABRIGHT, &
                              TABSTRING, TABUCSTRING, TABINTEGER, TABREAL
   use MemoryHelperModule, only: create_mem_path
-  use GenericUtilitiesModule, only: is_same, sim_message
+  use GenericUtilitiesModule, only: is_close, sim_message
   use SmoothingModule, only: sQuadraticSaturation, &
                              sQuadraticSaturationDerivative, &
                              sQuadratic0sp, &

src/Model/GroundWaterFlow/gwf3lak8.f90

@@ -27,7 +27,7 @@ module LakModule
   use BaseDisModule, only: DisBaseType
   use SimModule, only: count_errors, store_error, store_error_unit, &
                        deprecation_warning
-  use GenericUtilitiesModule, only: sim_message, is_same
+  use GenericUtilitiesModule, only: sim_message, is_close
   use BlockParserModule, only: BlockParserType
   use BaseDisModule, only: DisBaseType
   use SimVariablesModule, only: errmsg, warnmsg
@@ -861,7 +861,7 @@ subroutine lak_read_lake_connections(this)
                                    warnmsg, this%parser%GetUnit())
         case default
           read (keyword, *) rval
-          if (is_same(rval, DNODATA)) then
+          if (is_close(rval, DNODATA)) then
             is_lake_bed = .FALSE.
           else
             is_lake_bed = .TRUE.
@@ -1904,7 +1904,7 @@ subroutine lak_read_initial_attr(this)
           end if
           length = this%connlength(j)
         end if
-        if (is_same(this%bedleak(j), DNODATA)) then
+        if (is_close(this%bedleak(j), DNODATA)) then
           clb(j) = DNODATA
         else if (this%bedleak(j) > DZERO) then
           clb(j) = DONE / this%bedleak(j)
@@ -1916,7 +1916,7 @@ subroutine lak_read_initial_attr(this)
         else
           caq(j) = DZERO
         end if
-        if (is_same(this%bedleak(j), DNODATA)) then
+        if (is_close(this%bedleak(j), DNODATA)) then
           this%satcond(j) = area / caq(j)
         else if (clb(j) * caq(j) > DZERO) then
           this%satcond(j) = area / (clb(j) + caq(j))
@@ -1951,7 +1951,7 @@ subroutine lak_read_initial_attr(this)
           nn = this%cellid(j)
           area = this%warea(j)
           c1 = DZERO
-          if (is_same(clb(j), DNODATA)) then
+          if (is_close(clb(j), DNODATA)) then
             cbedleak = '     NONE     '
             cbedcond = '     NONE     '
           else if (clb(j) > DZERO) then
@@ -2767,7 +2767,7 @@ subroutine lak_calculate_evaporation(this, ilak, stage, avail, ev)
     call this%lak_calculate_sarea(ilak, stage, sa)
     ev = sa * this%evaporation(ilak)
     if (ev > avail) then
-      if (is_same(avail, DPREC)) then
+      if (is_close(avail, DPREC)) then
         ev = DZERO
       else
         ev = -avail

src/Solution/LinearMethods/ims8base.f90

@@ -9,7 +9,7 @@ MODULE IMSLinearBaseModule
   use KindModule, only: DP, I4B
   use ConstantsModule, only: LINELENGTH, IZERO, &
                              DZERO, DPREC, DEM6, DEM3, DHALF, DONE
-  use GenericUtilitiesModule, only: sim_message, is_same
+  use GenericUtilitiesModule, only: sim_message, is_close
   use BlockParserModule, only: BlockParserType
   use IMSReorderingModule, only: ims_odrv
   use ConvergenceSummaryModule
@@ -211,7 +211,7 @@ SUBROUTINE ims_base_cg(ICNVG, ITMAX, INNERIT, &
       IF (ICNVG .NE. 0) EXIT INNER
       !
       ! -- CHECK THAT CURRENT AND PREVIOUS rho ARE DIFFERENT
-      lsame = is_same(rho, rho0)
+      lsame = is_close(rho, rho0)
       IF (lsame) THEN
         EXIT INNER
       END IF
@@ -513,15 +513,15 @@ SUBROUTINE ims_base_bcgs(ICNVG, ITMAX, INNERIT, &
       !
       ! -- CHECK THAT CURRENT AND PREVIOUS rho, alpha, AND omega ARE
       !    DIFFERENT
-      lsame = is_same(rho, rho0)
+      lsame = is_close(rho, rho0)
       IF (lsame) THEN
         EXIT INNER
       END IF
-      lsame = is_same(alpha, alpha0)
+      lsame = is_close(alpha, alpha0)
       IF (lsame) THEN
         EXIT INNER
       END IF
-      lsame = is_same(omega, omega0)
+      lsame = is_close(omega, omega0)
       IF (lsame) THEN
         EXIT INNER
       END IF

src/Solution/NumericalSolution.f90

@@ -12,7 +12,7 @@ module NumericalSolutionModule
                              LENMEMPATH
   use MemoryHelperModule, only: create_mem_path
   use TableModule, only: TableType, table_cr
-  use GenericUtilitiesModule, only: is_same, sim_message, stop_with_error
+  use GenericUtilitiesModule, only: is_close, sim_message, stop_with_error
   use VersionModule, only: IDEVELOPMODE
   use BaseModelModule, only: BaseModelType
   use BaseExchangeModule, only: BaseExchangeType
@@ -2470,7 +2470,7 @@ subroutine sln_ls(this, kiter, kstp, kper, in_iter, iptc, ptcf)
           end if
         end if
       else
-        lsame = is_same(l2norm, this%l2norm0)
+        lsame = is_close(l2norm, this%l2norm0)
         if (lsame) then
           iptc = 0
         end if

src/Utilities/InputOutput.f90

@@ -11,7 +11,7 @@ module InputOutputModule
                              TABLEFT, TABCENTER, TABRIGHT,                     &
                              TABSTRING, TABUCSTRING, TABINTEGER, TABREAL,      &
                              DZERO
-  use GenericUtilitiesModule, only: is_same, sim_message
+  use GenericUtilitiesModule, only: is_close, sim_message
   private
   public :: GetUnit,                                                           &
             UPCASE, URWORD, ULSTLB, UBDSV4,                                    &

src/Utilities/MathUtil.f90

@@ -0,0 +1,49 @@
+!> @brief Generic math utilities
+module MathUtilitiesModule
+  use KindModule, only: DP, I4B, LGP
+  use ConstantsModule, only: DZERO, DPREC, DSAME
+  implicit none
+  private
+  public :: is_close
+
+contains
+
+  !> @brief Determines whether two double precision reals are close.
+  elemental function is_close(a, b, rel_tol, abs_tol) result(close)
+    ! -- return
+    logical(LGP) close
+    ! -- dummy
+    real(DP), intent(in) :: a, b
+    real(DP), intent(in), optional :: rel_tol, abs_tol
+    ! -- local
+    real(DP) :: loc_rel_tol, loc_abs_tol
+    !
+    ! -- evaluate optionals
+    if (present(rel_tol)) then
+      loc_rel_tol = rel_tol
+    else
+      loc_rel_tol = DSAME
+    end if
+    close = .false.
+    if (a == b) then
+      close = .true.
+    else
+      if (abs(b) > abs(a)) then
+        denom = b
+      else
+        denom = a
+        if (abs(denom) == DZERO) then
+          denom = DPREC
+        end if
+      end if
+      rdiff = abs((a - b) / denom)
+      if (rdiff <= epsloc) then
+        close = .true.
+      end if
+    end if
+    !
+    ! -- return
+    return
+  end function is_close
+
+end module MathUtilitiesModule

src/Utilities/TimeSeries/TimeArraySeries.f90

@@ -4,7 +4,7 @@ module TimeArraySeriesModule
   use BlockParserModule, only: BlockParserType
   use ConstantsModule, only: LINELENGTH, UNDEFINED, STEPWISE, LINEAR, &
                              LENTIMESERIESNAME, LENMODELNAME, DZERO, DONE
-  use GenericUtilitiesModule, only: is_same
+  use GenericUtilitiesModule, only: is_close
   use InputOutputModule, only: GetUnit, openfile
   use KindModule, only: DP, I4B
   use ListModule, only: ListType, ListNodeType
@@ -502,7 +502,7 @@ subroutine get_values_at_time(this, nvals, values, time)
           ierr = 1
         end if
       else
-        if (is_same(taEarlier%taTime, time)) then
+        if (is_close(taEarlier%taTime, time)) then
           values = taEarlier%taArray
         else
           ! -- Only earlier time is available, and it is not time of interest;
@@ -516,7 +516,7 @@ subroutine get_values_at_time(this, nvals, values, time)
       end if
     else
       if (associated(taLater)) then
-        if (is_same(taLater%taTime, time)) then
+        if (is_close(taLater%taTime, time)) then
           values = taLater%taArray
         else
           ! -- only later time is available, and it is not time of interest
@@ -756,7 +756,7 @@ subroutine get_latest_preceding_node(this, time, tslNode)
         if (associated(currNode%nextNode)) then
           obj => currNode%nextNode%GetItem()
           ta => CastAsTimeArrayType(obj)
-          if (ta%taTime < time .or. is_same(ta%taTime, time)) then
+          if (ta%taTime < time .or. is_close(ta%taTime, time)) then
             currNode => currNode%nextNode
           else
             exit