Merge pull request #488 from fortran-fans/add_close

[stdlib_math] add `is_close` routines.

Author: jvdp1

doc/specs/stdlib_math.md

@@ -320,25 +320,25 @@ program demo_logspace_rstart_cbase
 
 end program demo_logspace_rstart_cbase
 ```
-## `arange`
+### `arange`
 
-### Status
+#### Status
 
 Experimental
 
-### Class
+#### Class
 
 Pure function.
 
-### Description
+#### Description
 
 Creates a one-dimensional `array` of the `integer/real` type with fixed-spaced values of given spacing, within a given interval.
 
-### Syntax
+#### Syntax
 
 `result = [[stdlib_math(module):arange(interface)]](start [, end, step])`
 
-### Arguments
+#### Arguments
 
 All arguments should be the same type and kind.
 
@@ -354,18 +354,18 @@ The default `end` value is the inputted `start` value.
 This is an `intent(in)` and `optional` argument.   
 The default `step` value is `1`.
 
-#### Warning
+##### Warning
 If `step = 0`, the `step` argument will be corrected to `1/1.0` by the internal process of the `arange` function.   
 If `step < 0`, the `step` argument will be corrected to `abs(step)` by the internal process of the `arange` function. 
 
-### Return value
+#### Return value
 
 Returns a one-dimensional `array` of fixed-spaced values.
 
 For `integer` type arguments, the length of the result vector is `(end - start)/step + 1`.  
 For `real` type arguments, the length of the result vector is `floor((end - start)/step) + 1`.
 
-### Example
+#### Example
 
 ```fortran
 program demo_math_arange
@@ -388,3 +388,143 @@ program demo_math_arange
 
 end program demo_math_arange
 ```
+
+### `is_close`
+
+#### Description
+
+Returns a boolean scalar/array where two scalars/arrays are element-wise equal within a tolerance.
+
+```fortran
+!> For `real` type
+is_close(a, b, rel_tol, abs_tol) = abs(a - b) <= max(rel_tol*(abs(a), abs(b)), abs_tol)
+
+!> and for `complex` type
+is_close(a, b, rel_tol, abs_tol) = is_close(a%re, b%re, rel_tol, abs_tol) .and. &
+                                   is_close(a%im, b%im, rel_tol, abs_tol)
+```
+
+#### Syntax
+
+`bool = [[stdlib_math(module):is_close(interface)]] (a, b [, rel_tol, abs_tol, equal_nan])`
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Elemental function.
+
+#### Arguments
+
+Note: All `real/complex` arguments must have same `kind`.  
+If the value of `rel_tol/abs_tol` is negative (not recommended), 
+it will be corrected to `abs(rel_tol/abs_tol)` by the internal process of `is_close`.
+
+`a`: Shall be a `real/complex` scalar/array.
+This argument is `intent(in)`.
+
+`b`: Shall be a `real/complex` scalar/array.
+This argument is `intent(in)`.
+
+`rel_tol`: Shall be a `real` scalar/array.
+This argument is `intent(in)` and `optional`, which is `sqrt(epsilon(..))` by default.
+
+`abs_tol`: Shall be a `real` scalar/array.
+This argument is `intent(in)` and `optional`, which is `0.0` by default.
+
+`equal_nan`: Shall be a `logical` scalar/array.
+This argument is `intent(in)` and `optional`, which is `.false.` by default.
+Whether to compare `NaN` values as equal. If `.true.`, 
+`NaN` values in `a` will be considered equal to `NaN` values in `b`.
+
+#### Result value
+
+Returns a `logical` scalar/array.
+
+#### Example
+
+```fortran
+program demo_math_is_close
+
+    use stdlib_math,  only: is_close
+    use stdlib_error, only: check
+    real :: x(2) = [1, 2], y, NAN
+    
+    y   = -3
+    NAN = sqrt(y)
+    
+    print *, is_close(x,[real :: 1, 2.1])       !! [T, F]
+    print *, is_close(2.0, 2.1, abs_tol=0.1)    !! T
+    print *, NAN, is_close(2.0, NAN), is_close(2.0, NAN, equal_nan=.true.)   !! NAN, F, F
+    print *, is_close(NAN, NAN), is_close(NAN, NAN, equal_nan=.true.)        !! F, T
+        
+end program demo_math_is_close
+```
+
+### `all_close`
+
+#### Description
+
+Returns a boolean scalar where two arrays are element-wise equal within a tolerance.
+
+#### Syntax
+
+`bool = [[stdlib_math(module):all_close(interface)]] (a, b [, rel_tol, abs_tol, equal_nan])`
+
+#### Status
+
+Experimental.
+
+#### Class
+
+Pure function.
+
+#### Arguments
+
+Note: All `real/complex` arguments must have same `kind`.  
+If the value of `rel_tol/abs_tol` is negative (not recommended), 
+it will be corrected to `abs(rel_tol/abs_tol)` by the internal process of `all_close`.
+
+`a`: Shall be a `real/complex` array.
+This argument is `intent(in)`.
+
+`b`: Shall be a `real/complex` array.
+This argument is `intent(in)`.
+
+`rel_tol`: Shall be a `real` scalar.
+This argument is `intent(in)` and `optional`, which is `sqrt(epsilon(..))` by default.
+
+`abs_tol`: Shall be a `real` scalar.
+This argument is `intent(in)` and `optional`, which is `0.0` by default.
+
+`equal_nan`: Shall be a `logical` scalar.
+This argument is `intent(in)` and `optional`, which is `.false.` by default.
+Whether to compare `NaN` values as equal. If `.true.`, 
+`NaN` values in `a` will be considered equal to `NaN` values in `b`.
+
+#### Result value
+
+Returns a `logical` scalar.
+
+#### Example
+
+```fortran
+program demo_math_all_close
+
+    use stdlib_math,  only: all_close
+    use stdlib_error, only: check
+    real    :: x(2) = [1, 2], y, NAN
+    complex :: z(4, 4)
+    
+    y   = -3
+    NAN = sqrt(y)
+    z   = (1.0, 1.0)
+    
+    print *, all_close(z+cmplx(1.0e-11, 1.0e-11), z)     !! T
+    print *, NAN, all_close([NAN], [NAN]), all_close([NAN], [NAN], equal_nan=.true.) 
+                                                         !! NAN, F, T
+    
+end program demo_math_all_close
+```

src/CMakeLists.txt

@@ -40,6 +40,8 @@ set(fppFiles
     stdlib_math_linspace.fypp
     stdlib_math_logspace.fypp
     stdlib_math_arange.fypp
+    stdlib_math_is_close.fypp
+    stdlib_math_all_close.fypp
     stdlib_string_type.fypp
     stdlib_string_type_constructor.fypp
     stdlib_strings_to_string.fypp

src/Makefile.manual

@@ -35,8 +35,10 @@ SRCFYPP = \
         stdlib_stats_distribution_normal.fypp \
         stdlib_stats_var.fypp \
         stdlib_math.fypp \
-	stdlib_math_linspace.fypp \
-	stdlib_math_logspace.fypp \
+        stdlib_math_linspace.fypp \
+        stdlib_math_logspace.fypp \
+        stdlib_math_is_close.fypp \
+        stdlib_math_all_close.fypp \
         stdlib_string_type.fypp \
         stdlib_string_type_constructor.fypp \
         stdlib_strings.fypp \
@@ -81,7 +83,8 @@ $(SRCGEN): %.f90: %.fypp common.fypp
 # Fortran module dependencies
 f18estop.o: stdlib_error.o
 stdlib_ascii.o: stdlib_kinds.o
-stdlib_bitsets.o: stdlib_kinds.o
+stdlib_bitsets.o: stdlib_kinds.o \
+    stdlib_optval.o
 stdlib_bitsets_64.o: stdlib_bitsets.o
 stdlib_bitsets_large.o: stdlib_bitsets.o
 stdlib_error.o: stdlib_optval.o
@@ -104,7 +107,8 @@ stdlib_io_npy_save.o: \
     stdlib_error.o \
     stdlib_strings.o
 stdlib_linalg.o: \
-    stdlib_kinds.o
+    stdlib_kinds.o \
+    stdlib_optval.o
 stdlib_linalg_diag.o: \
     stdlib_linalg.o \
     stdlib_kinds.o
@@ -195,6 +199,11 @@ stdlib_math_logspace.o: \
 stdlib_math_arange.o: \
     stdlib_math.o
 stdlib_linalg_outer_product.o: stdlib_linalg.o
+stdlib_math_is_close.o: \
+    stdlib_math.o
+stdlib_math_all_close.o: \
+    stdlib_math.o \
+    stdlib_math_is_close.o
 stdlib_stringlist_type.o: stdlib_string_type.o \
     stdlib_math.o \
-	stdlib_optval.o
+    stdlib_optval.o

src/stdlib_math.fypp

@@ -14,7 +14,7 @@ module stdlib_math
     public :: EULERS_NUMBER_QP
 #:endif
     public :: DEFAULT_LINSPACE_LENGTH, DEFAULT_LOGSPACE_BASE, DEFAULT_LOGSPACE_LENGTH
-    public :: arange
+    public :: arange, is_close, all_close
 
     integer, parameter :: DEFAULT_LINSPACE_LENGTH = 100
     integer, parameter :: DEFAULT_LOGSPACE_LENGTH = 50
@@ -294,6 +294,39 @@ module stdlib_math
         #:endfor
     end interface arange
 
+    !> Version: experimental
+    !>
+    !> Returns a boolean scalar/array where two scalar/arrays are element-wise equal within a tolerance.
+    !> ([Specification](../page/specs/stdlib_math.html#is_close))
+    interface is_close
+        #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+        #:for k1, t1 in RC_KINDS_TYPES
+        elemental module logical function is_close_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
+            ${t1}$, intent(in) :: a, b
+            real(${k1}$), intent(in), optional :: rel_tol, abs_tol
+            logical, intent(in), optional :: equal_nan
+        end function is_close_${t1[0]}$${k1}$
+        #:endfor
+    end interface is_close
+
+    !> Version: experimental
+    !>
+    !> Returns a boolean scalar where two arrays are element-wise equal within a tolerance.
+    !> ([Specification](../page/specs/stdlib_math.html#all_close))
+    interface all_close
+        #:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+        #:set RANKS = range(1, MAXRANK + 1)
+        #:for k1, t1 in RC_KINDS_TYPES
+        #:for r1 in RANKS
+        logical pure module function all_close_${r1}$_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
+            ${t1}$, intent(in) :: a${ranksuffix(r1)}$, b${ranksuffix(r1)}$
+            real(${k1}$), intent(in), optional :: rel_tol, abs_tol
+            logical, intent(in), optional :: equal_nan
+        end function all_close_${r1}$_${t1[0]}$${k1}$
+        #:endfor
+        #:endfor
+    end interface all_close
+
 contains
 
     #:for k1, t1 in IR_KINDS_TYPES

src/stdlib_math_all_close.fypp

@@ -0,0 +1,25 @@
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+
+submodule (stdlib_math) stdlib_math_all_close
+
+    implicit none
+    
+contains
+
+    #:for k1, t1 in RC_KINDS_TYPES
+    #:for r1 in RANKS
+    logical pure module function all_close_${r1}$_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
+
+        ${t1}$, intent(in) :: a${ranksuffix(r1)}$, b${ranksuffix(r1)}$
+        real(${k1}$), intent(in), optional :: rel_tol, abs_tol
+        logical, intent(in), optional :: equal_nan
+
+        close = all(is_close(a, b, rel_tol, abs_tol, equal_nan))
+
+    end function all_close_${r1}$_${t1[0]}$${k1}$
+    #:endfor
+    #:endfor
+
+end submodule stdlib_math_all_close

src/stdlib_math_arange.fypp

@@ -44,7 +44,7 @@ contains
         step_  = optval(step, 1_${k1}$)
         step_  = sign(merge(step_, 1_${k1}$, step_ /= 0_${k1}$), end_ - start_)
 
-        allocate(result((end_ - start_)/step_ + 1))
+        allocate(result((end_ - start_)/step_ + 1_${k1}$))
 
         result = [(i, i=start_, end_, step_)]
 

src/stdlib_math_is_close.fypp

@@ -0,0 +1,47 @@
+#:include "common.fypp"
+
+submodule(stdlib_math) stdlib_math_is_close
+
+    use, intrinsic :: ieee_arithmetic, only: ieee_is_nan
+    implicit none
+    
+    #:for k1 in REAL_KINDS
+    real(${k1}$), parameter :: sqrt_eps_${k1}$ = sqrt(epsilon(1.0_${k1}$))
+    #:endfor
+    
+contains
+
+    #! Determines whether the values of `a` and `b` are close.
+
+    #:for k1, t1 in REAL_KINDS_TYPES
+    elemental module logical function is_close_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
+        ${t1}$, intent(in) :: a, b
+        real(${k1}$), intent(in), optional :: rel_tol, abs_tol
+        logical, intent(in), optional :: equal_nan
+        logical :: equal_nan_
+
+        equal_nan_ = optval(equal_nan, .false.)
+        
+        if (ieee_is_nan(a) .or. ieee_is_nan(b)) then
+            close = merge(.true., .false., equal_nan_ .and. ieee_is_nan(a) .and. ieee_is_nan(b))
+        else
+            close = abs(a - b) <= max( abs(optval(rel_tol, sqrt_eps_${k1}$)*max(abs(a), abs(b))), &
+                                       abs(optval(abs_tol, 0.0_${k1}$)) )
+        end if     
+
+    end function is_close_${t1[0]}$${k1}$
+    #:endfor
+
+    #:for k1, t1 in CMPLX_KINDS_TYPES
+    elemental module logical function is_close_${t1[0]}$${k1}$(a, b, rel_tol, abs_tol, equal_nan) result(close)
+        ${t1}$, intent(in) :: a, b
+        real(${k1}$), intent(in), optional :: rel_tol, abs_tol
+        logical, intent(in), optional :: equal_nan
+
+        close = is_close_r${k1}$(a%re, b%re, rel_tol, abs_tol, equal_nan) .and. &
+                is_close_r${k1}$(a%im, b%im, rel_tol, abs_tol, equal_nan)
+
+    end function is_close_${t1[0]}$${k1}$
+    #:endfor
+
+end submodule stdlib_math_is_close