refactor(MathUtil): remove implicits from root-finding fns (MODFLOW-U…

…SGS#1554)
geohackweek · Jan 17, 2024 · 314ecc5 · 314ecc5
1 parent 4cc6ac8
commit 314ecc5
Show file tree

Hide file tree

Showing 2 changed files with 45 additions and 17 deletions.
diff --git a/autotest/TestMathUtil.f90 b/autotest/TestMathUtil.f90
@@ -3,7 +3,7 @@ module TestMathUtil
   use ConstantsModule, only: DNODATA, DZERO
   use testdrive, only: check, error_type, new_unittest, test_failed, &
                        to_string, unittest_type
-  use MathUtilModule, only: is_close, mod_offset, &
+  use MathUtilModule, only: f1d, is_close, mod_offset, &
                             zeroch, zerotest, zeroin
   implicit none
   private
@@ -181,16 +181,19 @@ subroutine test_zeroch(error)
     type(error_type), allocatable, intent(out) :: error
     real(DP), parameter :: pi = 4 * atan(1.0_DP)
     real(DP) :: z
+    procedure(f1d), pointer :: f
 
-    z = zeroch(-1.0_DP, 1.0_DP, sine, 0.001_DP)
+    f => sine
+
+    z = zeroch(-1.0_DP, 1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, 0.0_DP, atol=1d-6), &
                'expected 0, got: '//to_string(z))
 
-    z = zeroch(-4.0_DP, -1.0_DP, sine, 0.001_DP)
+    z = zeroch(-4.0_DP, -1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, -pi, atol=1d-6), &
                'expected -pi, got: '//to_string(z))
 
-    z = zeroch(1.0_DP, 4.0_DP, sine, 0.001_DP)
+    z = zeroch(1.0_DP, 4.0_DP, f, 0.001_DP)
     call check(error, is_close(z, pi, atol=1d-6), &
                'expected pi, got: '//to_string(z))
   end subroutine test_zeroch
@@ -199,16 +202,19 @@ subroutine test_zeroin(error)
     type(error_type), allocatable, intent(out) :: error
     real(DP), parameter :: pi = 4 * atan(1.0_DP)
     real(DP) :: z
+    procedure(f1d), pointer :: f
+
+    f => sine
 
-    z = zeroin(-1.0_DP, 1.0_DP, sine, 0.001_DP)
+    z = zeroin(-1.0_DP, 1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, 0.0_DP, atol=1d-6), &
                'expected 0, got: '//to_string(z))
 
-    z = zeroin(-4.0_DP, -1.0_DP, sine, 0.001_DP)
+    z = zeroin(-4.0_DP, -1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, -pi, atol=1d-6), &
                'expected -pi, got: '//to_string(z))
 
-    z = zeroin(1.0_DP, 4.0_DP, sine, 0.001_DP)
+    z = zeroin(1.0_DP, 4.0_DP, f, 0.001_DP)
     call check(error, is_close(z, pi, atol=1d-6), &
                'expected pi, got: '//to_string(z))
   end subroutine test_zeroin
@@ -217,16 +223,19 @@ subroutine test_zerotest(error)
     type(error_type), allocatable, intent(out) :: error
     real(DP), parameter :: pi = 4 * atan(1.0_DP)
     real(DP) :: z
+    procedure(f1d), pointer :: f
+
+    f => sine
 
-    z = zerotest(-1.0_DP, 1.0_DP, sine, 0.001_DP)
+    z = zerotest(-1.0_DP, 1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, 0.0_DP, atol=1d-6), &
                'expected 0, got: '//to_string(z))
 
-    z = zerotest(-4.0_DP, -1.0_DP, sine, 0.001_DP)
+    z = zerotest(-4.0_DP, -1.0_DP, f, 0.001_DP)
     call check(error, is_close(z, -pi, atol=1d-6), &
                'expected -pi, got: '//to_string(z))
 
-    z = zerotest(1.0_DP, 4.0_DP, sine, 0.001_DP)
+    z = zerotest(1.0_DP, 4.0_DP, f, 0.001_DP)
     call check(error, is_close(z, pi, atol=1d-6), &
                'expected pi, got: '//to_string(z))
   end subroutine test_zerotest

diff --git a/src/Utilities/MathUtil.f90 b/src/Utilities/MathUtil.f90
@@ -7,12 +7,20 @@ module MathUtilModule
 
   implicit none
   private
-  public :: mod_offset, is_close, zeroch, zeroin, zerotest
+  public :: f1d, is_close, mod_offset, zeroch, zeroin, zerotest
 
   interface mod_offset
     module procedure :: mod_offset_int, mod_offset_dbl
   end interface mod_offset
 
+  interface
+    function f1d(x) result(fx)
+      import DP
+      real(DP), intent(in) :: x
+      real(DP) :: fx
+    end function
+  end interface
+
 contains
 
   !> @brief Check if a real value is approximately equal to another.
@@ -121,7 +129,11 @@ end function mod_offset_dbl
   !!
   !<
   function zeroch(x0, x1, f, epsa) result(z)
-    implicit double precision(a - h, o - z)
+    ! -- dummy
+    real(DP) :: x0, x1
+    procedure(f1d), pointer, intent(in) :: f
+    real(DP) :: epsa
+    real(DP) :: z
     ! -- local
     real(DP) :: epsm
     real(DP) :: a, b, c, t
@@ -131,7 +143,7 @@ function zeroch(x0, x1, f, epsa) result(z)
     real(DP) :: phi, philo, phihi
     real(DP) :: racb, rcab, rbca
     real(DP) :: tol, tl, tlc
-    real(DP) :: xm, xt
+    real(DP) :: xi, xm, xt
 
     epsm = epsilon(x0)
     b = x0
@@ -241,12 +253,15 @@ function zeroch(x0, x1, f, epsa) result(z)
   !! minimization without derivatives, prentice-hall, inc. (1973).
   !<
   function zeroin(ax, bx, f, tol) result(z)
-    implicit double precision(a - h, o - z)
+    ! -- dummy
+    real(DP) :: ax, bx
+    procedure(f1d), pointer, intent(in) :: f
+    real(DP) :: tol
+    real(DP) :: z
     ! -- local
     real(DP) :: eps
     real(DP) :: a, b, c, d, e, s, p, q
-    real(DP) :: fa, fb, fc
-    real(DP) :: tol1
+    real(DP) :: fa, fb, fc, r, tol1, xm
     logical(LGP) :: rs
 
     eps = epsilon(ax)
@@ -341,7 +356,11 @@ end function zeroin
 
   !> @brief Compute a zero of the function f(x) in the interval (x0, x1)
   function zerotest(x0, x1, f, epsa) result(z)
-    implicit double precision(a - h, o - z)
+    ! -- dummy
+    real(DP) :: x0, x1
+    procedure(f1d), pointer, intent(in) :: f
+    real(DP) :: epsa
+    real(DP) :: z
     ! -- local
     real(DP) :: epsm
     real(DP) :: ema, emb