diff --git a/doc/specs/stdlib_experimental_stats.md b/doc/specs/stdlib_experimental_stats.md
index cb955e5cc..f7bf46356 100644
--- a/doc/specs/stdlib_experimental_stats.md
+++ b/doc/specs/stdlib_experimental_stats.md
@@ -6,6 +6,54 @@ title: experimental_stats
 
 [TOC]
 
+## `corr` - Pearson correlation of array elements
+
+### Description
+
+Returns the Pearson correlation of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true`.
+
+The Pearson correlation between two rows (or columns), say `x` and `y`, of `array` is defined as:
+
+```
+ corr(x, y) = cov(x, y) / sqrt( var(x) * var(y))
+```
+
+### Syntax
+
+`result = [[stdlib_experimental_stats(module):corr(interface)]](array, dim [, mask])`
+
+### Arguments
+
+`array`: Shall be a rank-1 or a rank-2 array of type `integer`, `real`, or `complex`.
+
+`dim`: Shall be a scalar of type `integer` with a value in the range from 1 to `n`, where `n` is the rank of `array`.
+
+`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
+
+### Return value
+
+If `array` is of rank 1 and of type `real` or `complex`, the result is of type `real` and has the same kind as `array`.
+If `array` is of rank 2 and of type `real` or `complex`, the result is of the same type and kind as `array`.
+If `array` is of type `integer`, the result is of type `real(dp)`.
+
+If `array` is of rank 1 and of size larger than 1, a scalar equal to 1 is returned. Otherwise, IEEE `NaN` is returned.
+If `array` is of rank 2, a rank-2 array  with the corresponding correlations is returned.
+
+If `mask` is specified, the result is the Pearson correlation of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
+
+### Example
+
+```fortran
+program demo_corr
+    use stdlib_experimental_stats, only: corr
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    real :: y(1:2, 1:3) = reshape([ -1., 40., -3., 4., 10., 6. ], [ 2, 3])
+    print *, corr(x, 1)           !returns 1.
+    print *, corr(y, 2)           !returns reshape([ 1., -.32480, -.32480, 1. ], [ 2, 3])
+end program demo_corr
+```
+
 ## `cov` - covariance of array elements
 
 ### Description
@@ -233,4 +281,3 @@ program demo_var
     print *, var(y, 1, y > 3., corrected=.false.)    !returns [NaN, 0., 0.25]
 end program demo_var
 ```
-
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 51fd01d62..1d67eba77 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -7,6 +7,7 @@ set(fppFiles
     stdlib_experimental_linalg_diag.fypp
     stdlib_experimental_optval.fypp
     stdlib_experimental_stats.fypp
+    stdlib_experimental_stats_corr.fypp
     stdlib_experimental_stats_cov.fypp
     stdlib_experimental_stats_mean.fypp
     stdlib_experimental_stats_moment.fypp
diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp
index 6d3d12700..f3632a8f0 100644
--- a/src/stdlib_experimental_stats.fypp
+++ b/src/stdlib_experimental_stats.fypp
@@ -11,11 +11,102 @@ module stdlib_experimental_stats
   implicit none
   private
   ! Public API
-  public :: cov, mean, moment, var
+  public :: corr, cov, mean, moment, var
+
+
+  interface corr
+    !! Pearson correlation of array elements
+    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description))
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("corr",1, t1, k1)
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("corr",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("corr_mask",1, t1, k1)
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("corr_mask",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("corr",2, t1, k1)
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("corr",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                        , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("corr_mask",2, t1, k1)
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("corr_mask",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                            , merge(size(x, 1), size(x, 2), mask = 1<dim))
+      end function ${RName}$
+    #:endfor
+
+  end interface corr
+
 
   interface cov
     !! Covariance of array elements
-    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description))
+    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_1))
     #:for k1, t1 in RC_KINDS_TYPES
       #:set RName = rname("cov",1, t1, k1)
       module function ${RName}$(x, dim, mask, corrected) result(res)
@@ -112,7 +203,7 @@ module stdlib_experimental_stats
 
   interface mean
     !! Mean of array elements
-    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_1))
+    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_2))
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("mean_all",rank, t1, k1)
@@ -210,7 +301,7 @@ module stdlib_experimental_stats
 
   interface var
     !! Variance of array elements
-    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_3))
+    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_4))
 
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
@@ -317,7 +408,7 @@ module stdlib_experimental_stats
 
   interface moment
     !! Central moment of array elements
-    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_2))
+    !! ([Specification](../page/specs/stdlib_experimental_stats.html#description_3))
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
         #:set RName = rname("moment_all",rank, t1, k1)
diff --git a/src/stdlib_experimental_stats_corr.fypp b/src/stdlib_experimental_stats_corr.fypp
new file mode 100644
index 000000000..7924cf194
--- /dev/null
+++ b/src/stdlib_experimental_stats_corr.fypp
@@ -0,0 +1,311 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_experimental_stats) stdlib_experimental_stats_corr
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+  use stdlib_experimental_error, only: error_stop
+  use stdlib_experimental_linalg, only: diag
+  use stdlib_experimental_optval, only: optval
+  implicit none
+
+contains
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("corr",1, t1, k1)
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      real(${k1}$) :: res
+
+      if (.not.optval(mask, .true.) .or. size(x) < 2) then
+        res = ieee_value(1._${k1}$, ieee_quiet_nan)
+        return
+      end if
+
+      res = 1
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("corr",1, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      real(dp) :: res
+
+      if (.not.optval(mask, .true.) .or. size(x) < 2) then
+        res = ieee_value(1._dp, ieee_quiet_nan)
+        return
+      end if
+
+      res = 1
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("corr_mask",1, t1, k1)
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:)
+      real(${k1}$) :: res
+
+      if (count(mask) < 2) then
+        res = ieee_value(1._${k1}$, ieee_quiet_nan)
+        return
+      end if
+
+      res = 1
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("corr_mask",1, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:)
+      real(dp) :: res
+
+      if (count(mask) < 2) then
+        res = ieee_value(1._dp, ieee_quiet_nan)
+        return
+      end if
+
+      res = 1
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("corr",2, t1, k1)
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                          , merge(size(x, 1), size(x, 2), mask = 1<dim))
+
+      integer :: i, j
+      ${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
+      ${t1}$ :: center(size(x, 1),size(x, 2))
+
+      if (.not.optval(mask, .true.) .or. size(x) < 2) then
+        res = ieee_value(1._${k1}$, ieee_quiet_nan)
+        return
+      end if
+
+      mean_ = mean(x, dim)
+      select case(dim)
+        case(1)
+          do i = 1, size(x, 1)
+            center(i, :) = x(i, :) - mean_
+          end do
+          #:if t1[0] == 'r'
+            res = matmul( transpose(center), center)
+          #:else
+            res = matmul( transpose(conjg(center)), center)
+          #:endif
+        case(2)
+          do i = 1, size(x, 2)
+            center(:, i) = x(:, i) - mean_
+          end do
+          #:if t1[0] == 'r'
+            res = matmul( center, transpose(center))
+          #:else
+            res = matmul( center, transpose(conjg(center)))
+          #:endif
+        case default
+          call error_stop("ERROR (corr): wrong dimension")
+      end select
+
+      mean_ = 1 / sqrt(diag(res))
+      do i = 1, size(res, 1)
+        do j = 1, size(res, 2)
+          res(j, i) = res(j, i) * mean_(i) * mean_(j)
+        end do
+      end do
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("corr",2, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                      , merge(size(x, 1), size(x, 2), mask = 1<dim))
+
+      integer :: i, j
+      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
+      real(dp) :: center(size(x, 1),size(x, 2))
+
+      if (.not.optval(mask, .true.) .or. size(x) < 2) then
+        res = ieee_value(1._dp, ieee_quiet_nan)
+        return
+      end if
+
+      mean_ = mean(x, dim)
+      select case(dim)
+        case(1)
+          do i = 1, size(x, 1)
+            center(i, :) = real(x(i, :), dp) - mean_
+          end do
+          res = matmul( transpose(center), center)
+        case(2)
+          do i = 1, size(x, 2)
+            center(:, i) = real(x(:, i), dp) - mean_
+          end do
+          res = matmul( center, transpose(center))
+        case default
+          call error_stop("ERROR (corr): wrong dimension")
+      end select
+
+      mean_ = 1 / sqrt(diag(res))
+      do i = 1, size(res, 1)
+        do j = 1, size(res, 2)
+          res(j, i) = res(j, i) * mean_(i) * mean_(j)
+        end do
+      end do
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("corr_mask",2, t1, k1)
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:,:)
+      ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                          , merge(size(x, 1), size(x, 2), mask = 1<dim))
+
+      integer :: i, j
+      ${t1}$ :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      ${t1}$ :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+
+      select case(dim)
+        case(1)
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
+             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+
+              res(j, i) = dot_product( centerj_, centeri_)&
+               /sqrt(dot_product( centeri_, centeri_)*&
+                     dot_product( centerj_, centerj_))
+
+            end do
+          end do
+        case(2)
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
+             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask_)
+
+              res(j, i) = dot_product( centeri_, centerj_)&
+               /sqrt(dot_product( centeri_, centeri_)*&
+                     dot_product( centerj_, centerj_))
+            end do
+          end do
+        case default
+          call error_stop("ERROR (corr): wrong dimension")
+      end select
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("corr_mask",2, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:,:)
+      real(dp) :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                          , merge(size(x, 1), size(x, 2), mask = 1<dim))
+
+      integer :: i, j
+      real(dp) :: centeri_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      real(dp) :: centerj_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+      logical :: mask_(merge(size(x, 2), size(x, 1), mask = 1<dim))
+
+      select case(dim)
+        case(1)
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(:, i) .and. mask(:, j))
+             centeri_ = merge( x(:, i) - mean(x(:, i), mask = mask_),&
+                0._dp, mask_)
+             centerj_ = merge( x(:, j) - mean(x(:, j), mask = mask_),&
+                0._dp, mask_)
+
+             res(j, i) = dot_product( centerj_, centeri_)&
+                 /sqrt(dot_product( centeri_, centeri_)*&
+                       dot_product( centerj_, centerj_))
+
+            end do
+          end do
+        case(2)
+          do i = 1, size(res, 2)
+            do j = 1, size(res, 1)
+             mask_ = merge(.true., .false., mask(i, :) .and. mask(j, :))
+             centeri_ = merge( x(i, :) - mean(x(i, :), mask = mask_),&
+                0._dp, mask_)
+             centerj_ = merge( x(j, :) - mean(x(j, :), mask = mask_),&
+                0._dp, mask_)
+
+             res(j, i) = dot_product( centeri_, centerj_)&
+                 /sqrt(dot_product( centeri_, centeri_)*&
+                       dot_product( centerj_, centerj_))
+            end do
+          end do
+        case default
+          call error_stop("ERROR (corr): wrong dimension")
+      end select
+
+    end function ${RName}$
+  #:endfor
+
+
+end submodule
diff --git a/src/tests/stats/CMakeLists.txt b/src/tests/stats/CMakeLists.txt
index ce891bc4a..36ffc7aeb 100644
--- a/src/tests/stats/CMakeLists.txt
+++ b/src/tests/stats/CMakeLists.txt
@@ -1,3 +1,4 @@
+ADDTEST(corr)
 ADDTEST(cov)
 ADDTEST(mean)
 ADDTEST(moment)
diff --git a/src/tests/stats/test_corr.f90 b/src/tests/stats/test_corr.f90
new file mode 100644
index 000000000..a7ff0bd31
--- /dev/null
+++ b/src/tests/stats/test_corr.f90
@@ -0,0 +1,387 @@
+program test_corr
+    use stdlib_experimental_error, only: check
+    use stdlib_experimental_kinds, only: sp, dp, int32, int64
+    use stdlib_experimental_stats, only: corr
+    use,intrinsic :: ieee_arithmetic, only : ieee_is_nan
+    implicit none
+
+
+    real(sp), parameter :: sptol = 1000 * epsilon(1._sp)
+    real(dp), parameter :: dptol = 1000 * epsilon(1._dp)
+
+    real(dp) :: d1(5) = [1.0_dp, 2.0_dp, 3.0_dp, 4.0_dp, 5.0_dp]
+    real(dp) :: d(4, 3) = reshape([1._dp, 3._dp, 5._dp, 22._dp,&
+                                   3._dp, 4._dp, 6._dp, 20._dp,&
+                                   15._dp, 14._dp, 13._dp, 12._dp], [4, 3])
+
+    complex(dp) :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp, kind = dp),&
+                            cmplx(0.00000_dp, 1.44065_dp, kind = dp),&
+                            cmplx(1.26401_dp, 0.00000_dp, kind = dp),&
+                            cmplx(0.00000_dp, 0.88833_dp, kind = dp),&
+                            cmplx(1.14352_dp, 0.00000_dp, kind = dp)]
+    complex(dp) :: ds(2,3) = reshape([ cmplx(1._dp, 0._dp, kind = dp),&
+                            cmplx(0._dp, 2._dp, kind = dp),&
+                            cmplx(3._dp, 0._dp, kind = dp),&
+                            cmplx(0._dp, 4._dp, kind = dp),&
+                            cmplx(5._dp, 0._dp, kind = dp),&
+                            cmplx(0._dp, 6._dp, kind = dp)], [2, 3])
+
+
+    call test_sp(real(d1, sp), real(d, sp))
+
+    call test_dp(d1,d)
+
+    call test_int32(int(d1, int32) ,int(d, int32))
+
+    call test_int64(int(d1, int64) ,int(d, int64))
+
+    call test_csp(cmplx(cd1, kind = sp), cmplx(ds, kind = sp))
+
+    call test_cdp(cd1, ds)
+
+contains
+
+    subroutine test_sp(x, x2)
+        real(sp), intent(in) :: x(:)
+        real(sp), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, 1) - 1._sp) < sptol&
+            , 'sp check 1')
+        call check( ieee_is_nan(corr(x, 1, .false.))&
+            , 'sp check 2')
+        call check( ieee_is_nan(corr(x, 1, x == 1.)), 'sp check 3')
+        call check( abs(corr(x, 1, x < 5) - 1._sp) < sptol, 'sp check 4')
+        call check( ieee_is_nan(corr(x, 1, x < -999)), 'sp check 5')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = .false.)))&
+            , 'sp check 6')
+        call check( any(ieee_is_nan(corr(x2, 2, mask = .false.)))&
+            , 'sp check 7')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+            1._sp, 0.9994439103600_sp, -0.870544389237152_sp, 0.99944391036_sp,&
+            1._sp, -0.86261576629742_sp, -0.87054438923715_sp,  -0.862615766297428_sp,&
+            1._sp ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 8')
+        call check( all( abs( corr(x2, 2) - reshape([&
+             1._sp, 0.998742137866914_sp,  0.999846989517886_sp, -0.998337488459582_sp,&
+             0.998742137866914_sp, 1._sp, 0.999466429486246_sp, -0.99419162560192020_sp,&
+             0.999846989517886_sp, 0.999466429486246_sp, 1._sp, -0.99717646495273815_sp,&
+             -0.998337488459582_sp, -0.994191625601920_sp, -0.997176464952738_sp, 1._sp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 9')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = x2 < 10)))&
+            , 'sp check 10')
+        call check( all( abs( corr(x2, 1, mask = x2 < 22) - reshape([&
+              1._sp, 0.981980506061965_sp, -1._sp&
+              ,0.981980506061965_sp, 1._sp, -0.862615766297428_sp&
+              ,-1._sp, -0.862615766297428_sp, 1._sp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 11')
+        call check( all( abs( corr(x2, 2, mask = x2 < 22) - reshape([&
+            1._sp, 0.998742137866914_sp, 0.999846989517886_sp, -1._sp&
+            ,0.998742137866914_sp, 1._sp, 0.999466429486246_sp, -1._sp&
+            ,0.999846989517886_sp, 0.999466429486246_sp, 1._sp, -1._sp&
+            ,-1._sp, -1._sp, -1._sp, 1._sp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 12')
+
+
+        call check( all(abs(corr(x2, 1, mask = x2 < 1000) - corr(x2, 1))&
+            < sptol)&
+            , 'sp check 13')
+
+        call check( all(abs(corr(x2, 2, mask = x2 < 1000) - corr(x2, 2))&
+            < sptol)&
+            , 'sp check 14')
+
+    end subroutine test_sp
+
+    subroutine test_dp(x, x2)
+        real(dp), intent(in) :: x(:)
+        real(dp), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, 1) - 1._dp) < dptol&
+            , 'dp check 1')
+        call check( ieee_is_nan(corr(x, 1, .false.))&
+            , 'dp check 2')
+        call check( ieee_is_nan(corr(x, 1, x == 1.)), 'dp check 3')
+        call check( abs(corr(x, 1, x < 5) - 1._dp) < dptol, 'dp check 4')
+        call check( ieee_is_nan(corr(x, 1, x < -999)), 'dp check 5')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = .false.)))&
+            , 'dp check 6')
+        call check( any(ieee_is_nan(corr(x2, 2, mask = .false.)))&
+            , 'dp check 7')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+            1._dp, 0.9994439103600_dp, -0.870544389237152_dp, 0.99944391036_dp,&
+            1._dp, -0.86261576629742_dp, -0.87054438923715_dp,  -0.862615766297428_dp,&
+            1._dp ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 8')
+        call check( all( abs( corr(x2, 2) - reshape([&
+             1._dp, 0.998742137866914_dp,  0.999846989517886_dp, -0.998337488459582_dp,&
+             0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -0.99419162560192020_dp,&
+             0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -0.99717646495273815_dp,&
+             -0.998337488459582_dp, -0.994191625601920_dp, -0.997176464952738_dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 9')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = x2 < 10)))&
+            , 'dp check 10')
+        call check( all( abs( corr(x2, 1, mask = x2 < 22) - reshape([&
+              1._dp, 0.981980506061965_dp, -1._dp&
+              ,0.981980506061965_dp, 1._dp, -0.862615766297428_dp&
+              ,-1._dp, -0.862615766297428_dp, 1._dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 11')
+        call check( all( abs( corr(x2, 2, mask = x2 < 22) - reshape([&
+            1._dp, 0.998742137866914_dp, 0.999846989517886_dp, -1._dp&
+            ,0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -1._dp&
+            ,0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -1._dp&
+            ,-1._dp, -1._dp, -1._dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 12')
+
+        call check( all(abs(corr(x2, 1, mask = x2 < 1000) - corr(x2, 1))&
+            < dptol)&
+            , 'dp check 13')
+
+        call check( all(abs(corr(x2, 2, mask = x2 < 1000) - corr(x2, 2))&
+            < dptol)&
+            , 'dp check 14')
+
+    end subroutine test_dp
+
+    subroutine test_int32(x, x2)
+        integer(int32), intent(in) :: x(:)
+        integer(int32), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, 1) - 1._dp) < dptol&
+            , 'int32 check 1')
+        call check( ieee_is_nan(corr(x, 1, .false.))&
+            , 'int32 check 2')
+        call check( ieee_is_nan(corr(x, 1, x == 1.)), 'int32 check 3')
+        call check( abs(corr(x, 1, x < 5) - 1._dp) < dptol, 'int32 check 4')
+        call check( ieee_is_nan(corr(x, 1, x < -999)), 'int32 check 5')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = .false.)))&
+            , 'int32 check 6')
+        call check( any(ieee_is_nan(corr(x2, 2, mask = .false.)))&
+            , 'int32 check 7')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+            1._dp, 0.9994439103600_dp, -0.870544389237152_dp, 0.99944391036_dp,&
+            1._dp, -0.86261576629742_dp, -0.87054438923715_dp,  -0.862615766297428_dp,&
+            1._dp ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 8')
+        call check( all( abs( corr(x2, 2) - reshape([&
+             1._dp, 0.998742137866914_dp,  0.999846989517886_dp, -0.998337488459582_dp,&
+             0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -0.99419162560192020_dp,&
+             0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -0.99717646495273815_dp,&
+             -0.998337488459582_dp, -0.994191625601920_dp, -0.997176464952738_dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 9')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = x2 < 10)))&
+            , 'int32 check 10')
+        call check( all( abs( corr(x2, 1, mask = x2 < 22) - reshape([&
+              1._dp, 0.981980506061965_dp, -1._dp&
+              ,0.981980506061965_dp, 1._dp, -0.862615766297428_dp&
+              ,-1._dp, -0.862615766297428_dp, 1._dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 11')
+        call check( all( abs( corr(x2, 2, mask = x2 < 22) - reshape([&
+            1._dp, 0.998742137866914_dp, 0.999846989517886_dp, -1._dp&
+            ,0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -1._dp&
+            ,0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -1._dp&
+            ,-1._dp, -1._dp, -1._dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 12')
+
+        call check( all(abs(corr(x2, 1, mask = x2 < 1000) - corr(x2, 1))&
+            < dptol)&
+            , 'int32 check 13')
+
+        call check( all(abs(corr(x2, 2, mask = x2 < 1000) - corr(x2, 2))&
+            < dptol)&
+            , 'int32 check 14')
+
+    end subroutine test_int32
+
+    subroutine test_int64(x, x2)
+        integer(int64), intent(in) :: x(:)
+        integer(int64), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, 1) - 1._dp) < dptol&
+            , 'int64 check 1')
+        call check( ieee_is_nan(corr(x, 1, .false.))&
+            , 'int64 check 2')
+        call check( ieee_is_nan(corr(x, 1, x == 1.)), 'int64 check 3')
+        call check( abs(corr(x, 1, x < 5) - 1._dp) < dptol, 'int64 check 4')
+        call check( ieee_is_nan(corr(x, 1, x < -999)), 'int64 check 5')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = .false.)))&
+            , 'int64 check 6')
+        call check( any(ieee_is_nan(corr(x2, 2, mask = .false.)))&
+            , 'int64 check 7')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+            1._dp, 0.9994439103600_dp, -0.870544389237152_dp, 0.99944391036_dp,&
+            1._dp, -0.86261576629742_dp, -0.87054438923715_dp,  -0.862615766297428_dp,&
+            1._dp ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 8')
+        call check( all( abs( corr(x2, 2) - reshape([&
+             1._dp, 0.998742137866914_dp,  0.999846989517886_dp, -0.998337488459582_dp,&
+             0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -0.99419162560192020_dp,&
+             0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -0.99717646495273815_dp,&
+             -0.998337488459582_dp, -0.994191625601920_dp, -0.997176464952738_dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 9')
+
+        call check( any(ieee_is_nan(corr(x2, 1, mask = x2 < 10)))&
+            , 'int64 check 10')
+        call check( all( abs( corr(x2, 1, mask = x2 < 22) - reshape([&
+              1._dp, 0.981980506061965_dp, -1._dp&
+              ,0.981980506061965_dp, 1._dp, -0.862615766297428_dp&
+              ,-1._dp, -0.862615766297428_dp, 1._dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 11')
+        call check( all( abs( corr(x2, 2, mask = x2 < 22) - reshape([&
+            1._dp, 0.998742137866914_dp, 0.999846989517886_dp, -1._dp&
+            ,0.998742137866914_dp, 1._dp, 0.999466429486246_dp, -1._dp&
+            ,0.999846989517886_dp, 0.999466429486246_dp, 1._dp, -1._dp&
+            ,-1._dp, -1._dp, -1._dp, 1._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 12')
+
+        call check( all(abs(corr(x2, 1, mask = x2 < 1000) - corr(x2, 1))&
+            < dptol)&
+            , 'int64 check 13')
+
+        call check( all(abs(corr(x2, 2, mask = x2 < 1000) - corr(x2, 2))&
+            < dptol)&
+            , 'int64 check 14')
+
+    end subroutine test_int64
+
+    subroutine test_csp(x, x2)
+        complex(sp), intent(in) :: x(:)
+        complex(sp), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, dim=1) - 1._sp)  < sptol&
+            , 'csp check 1')
+        call check( abs(corr(x, 1, aimag(x) == 0) - 1._sp ) < sptol&
+            , 'csp check 2')
+
+        call check( ieee_is_nan(corr(x, 1, aimag(x) == -99 )) &
+            , 'csp check 3')
+
+        call check( ieee_is_nan(real(corr(x, 1, .false.)))&
+            , 'csp check 4')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+               (1._sp,0._sp), (0.983869910099907_sp,-0.178885438199983_sp),&
+               (0.973417168333576_sp,-0.229039333725547_sp),&
+               (0.983869910099907_sp,0.178885438199983_sp), (1._sp,0._sp),&
+               (0.998687663476588_sp,-0.051214751973158_sp),&
+               (0.973417168333575_sp,0.229039333725547_sp),&
+               (0.998687663476588_sp,0.0512147519731583_sp), (1._sp,0._sp) ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'csp check 6')
+        call check( all( abs( corr(x2, 2) - reshape([&
+            (1._sp,0._sp), (0._sp,1._sp),&
+            (0._sp,-1._sp), (1._sp,0._sp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'csp check 7')
+
+        call check( all( abs( corr(x2, 2, mask = aimag(x2) < 6) - reshape([&
+             (1._sp,0._sp), (0._sp,1._sp)&
+             ,(0._sp,-1._sp), (1._sp,0._sp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'csp check 8')
+
+        call check( all(abs(corr(x2, 1, mask = aimag(x2) < 1000) - corr(x2, 1))&
+            < sptol)&
+            , 'csp check 9')
+
+        call check( all(abs(corr(x2, 2, mask = aimag(x2) < 1000) - corr(x2, 2))&
+            < sptol)&
+            , 'csp check 10')
+
+    end subroutine test_csp
+
+    subroutine test_cdp(x, x2)
+        complex(dp), intent(in) :: x(:)
+        complex(dp), intent(in) :: x2(:, :)
+
+        call check( abs(corr(x, dim=1) - 1._dp)  < dptol&
+            , 'cdp check 1')
+        call check( abs(corr(x, 1, aimag(x) == 0) - 1._dp ) < dptol&
+            , 'cdp check 2')
+
+        call check( ieee_is_nan(corr(x, 1, aimag(x) == -99 )) &
+            , 'cdp check 3')
+
+        call check( ieee_is_nan(real(corr(x, 1, .false.)))&
+            , 'cdp check 4')
+
+        call check( all( abs( corr(x2, 1) - reshape([&
+               (1._dp,0._dp), (0.983869910099907_dp,-0.178885438199983_dp),&
+               (0.973417168333576_dp,-0.229039333725547_dp),&
+               (0.983869910099907_dp,0.178885438199983_dp), (1._dp,0._dp),&
+               (0.998687663476588_dp,-0.051214751973158_dp),&
+               (0.973417168333575_dp,0.229039333725547_dp),&
+               (0.998687663476588_dp,0.0512147519731583_dp), (1._dp,0._dp) ]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'cdp check 6')
+        call check( all( abs( corr(x2, 2) - reshape([&
+            (1._dp,0._dp), (0._dp,1._dp),&
+            (0._dp,-1._dp), (1._dp,0._dp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 7')
+
+        call check( all( abs( corr(x2, 2, mask = aimag(x2) < 6) - reshape([&
+             (1._dp,0._dp), (0._dp,1._dp)&
+             ,(0._dp,-1._dp), (1._dp,0._dp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 8')
+
+        call check( all(abs(corr(x2, 1, mask = aimag(x2) < 1000) - corr(x2, 1))&
+            < sptol)&
+            , 'csp check 9')
+
+        call check( all(abs(corr(x2, 2, mask = aimag(x2) < 1000) - corr(x2, 2))&
+            < sptol)&
+            , 'csp check 10')
+
+    end subroutine test_cdp
+
+end program test_corr