diff --git a/doc/specs/stdlib_math.md b/doc/specs/stdlib_math.md
index e604a06b1..2ce8bfce6 100644
--- a/doc/specs/stdlib_math.md
+++ b/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
+```
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b0c969004..239b10833 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -39,6 +39,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
diff --git a/src/Makefile.manual b/src/Makefile.manual
index befa540ae..3c6105583 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -34,8 +34,10 @@ SRCFYPP = \
         stdlib_stats_distribution_uniform.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 \
@@ -80,7 +82,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
@@ -103,7 +106,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
@@ -189,6 +193,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
diff --git a/src/stdlib_math.fypp b/src/stdlib_math.fypp
index 89f412ff0..2db5d2a56 100644
--- a/src/stdlib_math.fypp
+++ b/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
diff --git a/src/stdlib_math_all_close.fypp b/src/stdlib_math_all_close.fypp
new file mode 100644
index 000000000..788587e62
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/src/stdlib_math_arange.fypp b/src/stdlib_math_arange.fypp
index 075d76f7a..eb3be3202 100644
--- a/src/stdlib_math_arange.fypp
+++ b/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_)]
 
diff --git a/src/stdlib_math_is_close.fypp b/src/stdlib_math_is_close.fypp
new file mode 100644
index 000000000..b97134802
--- /dev/null
+++ b/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
diff --git a/src/tests/math/test_stdlib_math.fypp b/src/tests/math/test_stdlib_math.fypp
index 651fea610..72c3adf5b 100644
--- a/src/tests/math/test_stdlib_math.fypp
+++ b/src/tests/math/test_stdlib_math.fypp
@@ -4,7 +4,7 @@
 
 module test_stdlib_math
     use testdrive, only : new_unittest, unittest_type, error_type, check, skip_test
-    use stdlib_math, only: clip
+    use stdlib_math, only: clip, is_close, all_close
     use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, xdp, qp
     implicit none
 
@@ -32,6 +32,14 @@ contains
             new_unittest("clip-real-double-bounds", test_clip_rdp_bounds), &
             new_unittest("clip-real-quad", test_clip_rqp), &
             new_unittest("clip-real-quad-bounds", test_clip_rqp_bounds) &
+            
+            !> Tests for `is_close` and `all_close`
+            #:for k1 in REAL_KINDS
+            , new_unittest("is_close-real-${k1}$", test_is_close_real_${k1}$) &
+            , new_unittest("is_close-cmplx-${k1}$", test_is_close_cmplx_${k1}$) &
+            , new_unittest("all_close-real-${k1}$", test_all_close_real_${k1}$) &
+            , new_unittest("all_close-cmplx-${k1}$", test_all_close_cmplx_${k1}$) &
+            #:endfor
             ]
 
     end subroutine collect_stdlib_math
@@ -203,6 +211,89 @@ contains
 #:endif
 
     end subroutine test_clip_rqp_bounds
+    
+    #:for k1 in REAL_KINDS
+    subroutine test_is_close_real_${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        real(${k1}$) :: x, NAN
+        x = -3; NAN = sqrt(x)
+        
+        call check(error, is_close(2.5_${k1}$, 2.5_${k1}$), .true.)
+        if (allocated(error)) return
+        call check(error, is_close(0.0_${k1}$, -0.0_${k1}$), .true.)
+        if (allocated(error)) return
+        call check(error, is_close(2.5_${k1}$, 1.2_${k1}$), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(NAN, NAN), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(NAN, 0.0_${k1}$), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(NAN, NAN, equal_nan=.true.), .true.)
+    
+    end subroutine test_is_close_real_${k1}$
+    
+    subroutine test_is_close_cmplx_${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        real(${k1}$) :: x, NAN
+        x = -3; NAN = sqrt(x)
+        
+        call check(error, is_close((2.5_${k1}$, 1.5_${k1}$), (2.5_${k1}$, 1.5_${k1}$)), .true.)
+        if (allocated(error)) return
+        call check(error, is_close((2.5_${k1}$, 1.2_${k1}$), (2.5_${k1}$, 1.5_${k1}$)), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(cmplx(NAN, NAN, ${k1}$), cmplx(NAN, NAN, ${k1}$)), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(cmplx(NAN, NAN, ${k1}$), cmplx(NAN, 0.0_${k1}$, ${k1}$)), .false.)
+        if (allocated(error)) return
+        call check(error, is_close(cmplx(NAN, NAN, ${k1}$), cmplx(NAN, NAN, ${k1}$), equal_nan=.true.), .true.)
+        if (allocated(error)) return
+        call check(error, is_close(cmplx(NAN, 1.2_${k1}$, ${k1}$), cmplx(NAN, 1.2_${k1}$, ${k1}$), equal_nan=.true.), .true.)
+        
+    end subroutine test_is_close_cmplx_${k1}$
+    
+    subroutine test_all_close_real_${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        real(${k1}$) :: x(2, 2), eps, NAN
+        x = 1; eps = -3; NAN = sqrt(eps)
+        
+        eps = sqrt(epsilon(1.0_${k1}$))
+        
+        call check(error, all_close(x, x), .true.)
+        if (allocated(error)) return
+        call check(error, all_close(x + x*eps + 1.0e-6, x), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + NAN, x), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + NAN, x, equal_nan=.true.), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + NAN, x + NAN), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + NAN, x + NAN, equal_nan=.true.), .true.)
+        
+    end subroutine test_all_close_real_${k1}$
+    
+    subroutine test_all_close_cmplx_${k1}$(error)
+        type(error_type), allocatable, intent(out) :: error
+        real(${k1}$) :: eps, NAN
+        complex(${k1}$) :: x(2, 2)
+        x = (1, 1); eps = -3; NAN = sqrt(eps)
+        
+        eps = sqrt(epsilon(1.0_${k1}$))
+        
+        call check(error, all_close(x, x), .true.)
+        if (allocated(error)) return
+        call check(error, all_close(x + x*eps + 1.0e-6, x), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + cmplx(NAN, NAN, ${k1}$), x), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + cmplx(NAN, NAN, ${k1}$), x, equal_nan=.true.), .false.)
+        if (allocated(error)) return
+        call check(error, all_close(x + cmplx(NAN, NAN, ${k1}$), x + cmplx(NAN, NAN, ${k1}$), equal_nan=.true.), .true.)
+        if (allocated(error)) return
+        call check(error, all_close(x + cmplx(NAN, NAN, ${k1}$), x + cmplx(NAN, NAN, ${k1}$)), .false.)
+        
+    end subroutine test_all_close_cmplx_${k1}$
+    #:endfor
 
 end module test_stdlib_math