From 285dab445192e4c78233b3a3f795a1aebca6b95d Mon Sep 17 00:00:00 2001
From: Martin Diehl <mail@martin-diehl.net>
Date: Sun, 22 Sep 2019 23:53:56 -0700
Subject: [PATCH] self consistent solutions

---
 src/math.f90 |  1 -
 src/prec.f90 | 62 +++++++++++++++++++++++++++-------------------------
 2 files changed, 32 insertions(+), 31 deletions(-)

diff --git a/src/math.f90 b/src/math.f90
index 174f8c6a7..9ffd0c325 100644
--- a/src/math.f90
+++ b/src/math.f90
@@ -1411,7 +1411,6 @@ subroutine unitTest
   real(pReal), dimension(9,9) :: t99,t99_2
   logical                     :: e
   
-integer :: i
   if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - &
               math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) &
     call IO_error(401,ext_msg='math_expand [1,2,3] by [1,2,3,0] => [1,2,2,3,3,3]')
diff --git a/src/prec.f90 b/src/prec.f90
index 2a1c2f14a..434ce32e7 100644
--- a/src/prec.f90
+++ b/src/prec.f90
@@ -10,6 +10,7 @@ module prec
   use, intrinsic :: IEEE_arithmetic
 
   implicit none
+  public
   ! https://software.intel.com/en-us/blogs/2017/03/27/doctor-fortran-in-it-takes-all-kinds
 #ifdef Abaqus
   integer,     parameter, public :: pReal      = selected_real_kind(15,307)                         !< number with 15 significant digits, up to 1e+-307 (typically 64 bit)
@@ -83,14 +84,8 @@ module prec
   real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal)                             !< minimum positive number such that 1.0 + EPSILON /= 1.0.
   real(pReal), private, parameter :: PREAL_MIN     = tiny(0.0_pReal)                                !< smallest normalized floating point number
 
-  public :: &
-    prec_init, &
-    dEq, &
-    dEq0, &
-    cEq, &
-    dNeq, &
-    dNeq0, &
-    cNeq
+  private :: &
+    unitTest
 
 contains
 
@@ -100,11 +95,6 @@ contains
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
 
-  integer, allocatable, dimension(:) :: realloc_lhs_test
-
-  external :: &
-    quit
-
   write(6,'(/,a)') ' <<<+-  prec init  -+>>>'
 
   write(6,'(a,i3)')    ' Size of integer in bit: ',bit_size(0)
@@ -114,8 +104,7 @@ subroutine prec_init
   write(6,'(a,e10.3)') '   Minimum value:        ',tiny(0.0_pReal)
   write(6,'(a,i3)')    '   Decimal precision:    ',precision(0.0_pReal)
 
-  realloc_lhs_test = [1,2]
-  if (realloc_lhs_test(2)/=2) call quit(9000)
+  call unitTest
 
 end subroutine prec_init
 
@@ -153,16 +142,13 @@ logical elemental pure function dNeq(a,b,tol)
 
   real(pReal), intent(in)           :: a,b
   real(pReal), intent(in), optional :: tol
-  real(pReal)                       :: eps
 
   if (present(tol)) then
-    eps = tol
+    dNeq = .not. dEq(a,b,tol)
   else
-    eps = PREAL_EPSILON * maxval(abs([a,b]))
+    dNeq = .not. dEq(a,b)
   endif
 
-  dNeq = merge(.False.,.True.,abs(a-b) <= eps)
-
 end function dNeq
 
 
@@ -199,16 +185,13 @@ logical elemental pure function dNeq0(a,tol)
 
   real(pReal), intent(in)           :: a
   real(pReal), intent(in), optional :: tol
-  real(pReal)                       :: eps
 
   if (present(tol)) then
-    eps = tol
+    dNeq0 = .not. dEq0(a,tol)
   else
-    eps = PREAL_MIN * 10.0_pReal
+    dNeq0 = .not. dEq0(a)
   endif
 
-  dNeq0 = merge(.False.,.True.,abs(a) <= eps)
-
 end function dNeq0
 
 
@@ -247,16 +230,35 @@ logical elemental pure function cNeq(a,b,tol)
 
   complex(pReal), intent(in)           :: a,b
   real(pReal),    intent(in), optional :: tol
-  real(pReal)                          :: eps
 
   if (present(tol)) then
-    eps = tol
+    cNeq = .not. cEq(a,b,tol)
   else
-    eps = PREAL_EPSILON * maxval(abs([a,b]))
+    cNeq = .not. cEq(a,b)
   endif
 
-  cNeq = merge(.False.,.True.,abs(a-b) <= eps)
-
 end function cNeq
 
+
+!--------------------------------------------------------------------------------------------------
+!> @brief check correctness of (some) prec functions
+!--------------------------------------------------------------------------------------------------
+subroutine unitTest
+
+  integer, allocatable, dimension(:) :: realloc_lhs_test
+  real(pReal), dimension(2) :: r
+  external :: &
+    quit
+  
+  call random_number(r)
+  r = r/minval(r)
+  if(.not. all(dEq(r,r+PREAL_EPSILON)))    call quit(9000)
+  if(dEq(r(1),r(2)) .and. dNeq(r(1),r(2))) call quit(9000)
+  if(.not. all(dEq0(r-r+PREAL_MIN)))       call quit(9000)
+  
+  realloc_lhs_test = [1,2]
+  if (any(realloc_lhs_test/=[1,2])) call quit(9000)
+  
+end subroutine unitTest
+
 end module prec