From 279a192041d7d3c6d2a3e82b50ee0203291115cb Mon Sep 17 00:00:00 2001
From: "Vandenplas, Jeremie" <jeremie.vandenplas@gmail.com>
Date: Sun, 12 Apr 2020 11:03:04 +0200
Subject: [PATCH 1/3] Init of covariance

Squashed commit of the following:

commit c7d6c1d501c151358c0fe8dc6bcabfe0424b95f9
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Sun Apr 12 10:18:15 2020 +0200

    covariance: addition of tests

commit 9cd7e63c50433d4c7901e64edfdfd9ec3a2df424
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Sun Apr 12 00:12:41 2020 +0200

    covariance:comment

commit 958b59be4f68a4ca74f6f048c52c1d0cbe45b209
Author: Vandenplas, Jeremie <jeremie.vandenplas@wur.nl>
Date:   Sat Apr 11 23:46:40 2020 +0200

    covariance: correct test

commit fdf4e54bc7a0958add6e9f99507d240e31c52ddc
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Sat Apr 11 22:57:23 2020 +0200

    covariance: update spec

commit c3c277f3aea79f4671994f5eb36ecff193638329
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Sat Apr 11 22:35:11 2020 +0200

    covariance: support of cov for vectors

commit c60ce03a7d74e4a9916e5eb21b1caf07a193bf86
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Sat Apr 11 21:46:41 2020 +0200

    covariance: real + complex+integer

commit 9e3eba041c07268c4782779744f4d434f20578d2
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Fri Apr 10 23:16:48 2020 +0200

    covariance: addition of complex

commit 340309e2ee0b8ddc29135a41876232ea5ee5f4e5
Author: Vandenplas, Jeremie <jeremie.vandenplas@gmail.com>
Date:   Fri Apr 10 22:32:17 2020 +0200

    covariance: init
---
 src/CMakeLists.txt                     |   1 +
 src/stdlib_experimental_stats.fypp     |  98 ++++-
 src/stdlib_experimental_stats.md       |  57 +++
 src/stdlib_experimental_stats_cov.fypp | 290 +++++++++++++
 src/tests/stats/CMakeLists.txt         |   1 +
 src/tests/stats/test_cov.f90           | 566 +++++++++++++++++++++++++
 6 files changed, 1012 insertions(+), 1 deletion(-)
 create mode 100644 src/stdlib_experimental_stats_cov.fypp
 create mode 100644 src/tests/stats/test_cov.f90

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 043e6932d..7e54cbeea 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -5,6 +5,7 @@ set(fppFiles
     stdlib_experimental_io.fypp
     stdlib_experimental_optval.fypp
     stdlib_experimental_stats.fypp
+    stdlib_experimental_stats_cov.fypp
     stdlib_experimental_stats_mean.fypp
     stdlib_experimental_stats_moment.fypp
     stdlib_experimental_stats_var.fypp
diff --git a/src/stdlib_experimental_stats.fypp b/src/stdlib_experimental_stats.fypp
index 4f3d7e6bb..e5cb89a2a 100644
--- a/src/stdlib_experimental_stats.fypp
+++ b/src/stdlib_experimental_stats.fypp
@@ -8,7 +8,102 @@ module stdlib_experimental_stats
   implicit none
   private
   ! Public API
-  public :: mean, moment, var
+  public :: cov, mean, moment, var
+
+  interface cov
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov",1, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov_mask",1, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        logical, intent(in), optional :: corrected
+        real(${k1}$) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in INT_KINDS_TYPES
+      #:set RName = rname("cov_mask",1, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:)
+        logical, intent(in), optional :: corrected
+        real(dp) :: res
+      end function ${RName}$
+    #:endfor
+
+    #:for k1, t1 in RC_KINDS_TYPES
+      #:set RName = rname("cov",2, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        ${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("cov",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in), optional :: mask
+        logical, intent(in), optional :: corrected
+        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("cov_mask",2, t1, k1)
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        logical, intent(in), optional :: corrected
+        ${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("cov_mask",2, t1, k1, 'dp')
+      module function ${RName}$(x, dim, mask, corrected) result(res)
+        ${t1}$, intent(in) :: x(:, :)
+        integer, intent(in) :: dim
+        logical, intent(in) :: mask(:,:)
+        logical, intent(in), optional :: corrected
+        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 cov
+
 
   interface mean
     #:for k1, t1 in RC_KINDS_TYPES
@@ -210,6 +305,7 @@ module stdlib_experimental_stats
 
   end interface var
 
+
   interface moment
     #:for k1, t1 in RC_KINDS_TYPES
       #:for rank in RANKS
diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md
index 2f6c77264..671351715 100644
--- a/src/stdlib_experimental_stats.md
+++ b/src/stdlib_experimental_stats.md
@@ -1,9 +1,66 @@
 # Descriptive statistics
 
+* [`cov` - covariance of array elements](#cov---covariance-of-array-elements)
 * [`mean` - mean of array elements](#mean---mean-of-array-elements)
 * [`moment` - central moments of array elements](#moment---central-moments-of-array-elements)
 * [`var` - variance of array elements](#var---variance-of-array-elements)
 
+
+## `cov` - covariance of array elements
+
+### Description
+
+Returns the covariance of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true` (if `mask` is provided).
+
+Per default, the covariance is defined as:
+
+```
+ cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)) )
+```
+
+where n is the number of elements.
+
+The scaling can be changed with the logical argument `corrected`. If `corrected` is `.false.`, then the sum is scaled with `n`, otherwise with `n-1`.
+
+
+### Syntax
+
+`result = cov(array, dim [, mask [, corrected]])`
+
+### Arguments
+
+`array`: Shall be a 1-rank or a 2-rank 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 by a scalar or an array of the same shape as `array`.
+
+`corrected` (optional): Shall be a scalar of type `logical`. If `corrected` is `.true.` (default value), the sum is scaled with `n-1`. If `corrected` is `.false.`, then the sum is scaled with `n`.
+
+### Return value
+
+If `array` is of rank 1 and of type `real` or `complex`, the result is of the same type `real` corresponding to the type of `array`.
+If `array` is of rank 2 and of type `real` or `complex`, the result is of the same type as `array`.
+If `array` is of type `integer`, the result is of type `real(dp)`.
+
+If `array` is of rank 1, a scalar with the covariance (that is the variance) of all elements in `array` is returned.
+If `array` is of rank 2, a 2-rank array is returned.
+
+If `mask` is specified, the result is the covariance 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_cov
+    use stdlib_experimental_stats, only: cov
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    real :: y(1:2, 1:3) = reshape([ 1., 2., 3., 4., 5., 6. ], [ 2, 3])
+    print *, cov(x, 1)                          !returns 3.5
+    print *, cov(x, 1, corrected = .false.)     !returns 2.9167
+    print *, cov(y, 1)                          !returns a square matrix of size 3 with all elements equal to 0.5
+end program demo_cov
+```
+
 ## `mean` - mean of array elements
 
 ### Description
diff --git a/src/stdlib_experimental_stats_cov.fypp b/src/stdlib_experimental_stats_cov.fypp
new file mode 100644
index 000000000..e19c48489
--- /dev/null
+++ b/src/stdlib_experimental_stats_cov.fypp
@@ -0,0 +1,290 @@
+#:include "common.fypp"
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_experimental_stats) stdlib_experimental_stats_cov
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+  use stdlib_experimental_error, only: error_stop
+  use stdlib_experimental_optval, only: optval
+  implicit none
+
+contains
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("cov",1, t1, k1)
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      logical, intent(in), optional :: corrected
+      real(${k1}$) :: res
+
+      if (.not.optval(mask, .true.)) then
+        res = ieee_value(1._${k1}$, ieee_quiet_nan)
+        return
+      end if
+
+      res = var(x, dim, corrected = optval(corrected, .true.))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("cov",1, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      logical, intent(in), optional :: corrected
+      real(dp) :: res
+
+      if (.not.optval(mask, .true.)) then
+        res = ieee_value(1._dp, ieee_quiet_nan)
+        return
+      end if
+
+      res = var(x, dim, corrected = optval(corrected, .true.))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("cov_mask",1, t1, k1)
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:)
+      logical, intent(in), optional :: corrected
+      real(${k1}$) :: res
+
+      res = var(x, dim, mask, corrected = optval(corrected, .true.))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("cov_mask",1, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:)
+      logical, intent(in), optional :: corrected
+      real(dp) :: res
+
+      res = var(x, dim, mask, corrected = optval(corrected, .true.))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("cov",2, t1, k1)
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      logical, intent(in), optional :: corrected
+      ${t1}$ :: res(merge(size(x, 1), size(x, 2), mask = 1<dim)&
+                          , merge(size(x, 1), size(x, 2), mask = 1<dim))
+
+      integer :: i
+      ${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.)) 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 (cov): wrong dimension")
+      end select
+      res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("cov",2, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in), optional :: mask
+      logical, intent(in), optional :: corrected
+      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
+      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.)) 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 (cov): wrong dimension")
+      end select
+      res = res / (size(x, dim) - merge(1, 0, optval(corrected, .true.)))
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:set RName = rname("cov_mask",2, t1, k1)
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:,:)
+      logical, intent(in), optional :: corrected
+      ${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, n
+      ${t1}$ :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
+      ${t1}$ :: center(size(x, 1),size(x, 2))
+
+      mean_ = mean(x, dim, mask = mask)
+      select case(dim)
+        case(1)
+          do i = 1, size(x, 1)
+            center(i, :) = merge( x(i, :) - mean_,&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask(i, :))
+          end do
+          #:if t1[0] == 'r'
+            res = matmul( transpose(center), center)
+          #:else
+            res = matmul( transpose(conjg(center)), center)
+          #:endif
+          do j = 1, size(res, 2)
+            do i = 1, size(res, 1)
+              n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
+              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                            optval(corrected, .true.) .and. n > 0))
+            end do
+          end do
+        case(2)
+          do i = 1, size(x, 2)
+            center(:, i) = merge( x(:, i) - mean_,&
+              #:if t1[0] == 'r'
+                0._${k1}$,&
+              #:else
+                cmplx(0,0,kind=${k1}$),&
+              #:endif
+                mask(:, i))
+          end do
+          #:if t1[0] == 'r'
+            res = matmul( center, transpose(center))
+          #:else
+            res = matmul( center, transpose(conjg(center)))
+          #:endif
+          do j = 1, size(res, 2)
+            do i = 1, size(res, 1)
+              n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
+              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                            optval(corrected, .true.) .and. n > 0))
+            end do
+          end do
+        case default
+          call error_stop("ERROR (cov): wrong dimension")
+      end select
+
+    end function ${RName}$
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:set RName = rname("cov_mask",2, t1, k1, 'dp')
+    module function ${RName}$(x, dim, mask, corrected) result(res)
+      ${t1}$, intent(in) :: x(:, :)
+      integer, intent(in) :: dim
+      logical, intent(in) :: mask(:,:)
+      logical, intent(in), optional :: corrected
+      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, n
+      real(dp) :: mean_(merge(size(x, 1), size(x, 2), mask = 1<dim))
+      real(dp) :: center(size(x, 1),size(x, 2))
+
+      mean_ = mean(x, dim, mask = mask)
+      select case(dim)
+        case(1)
+          do i = 1, size(x, 1)
+            center(i, :) = merge( x(i, :) - mean_,&
+                0._dp,&
+                mask(i, :))
+          end do
+          res = matmul( transpose(center), center)
+          do j = 1, size(res, 2)
+            do i = 1, size(res, 1)
+              n = count(merge(.true., .false., mask(:, i) .and. mask(:, j)))
+              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                            optval(corrected, .true.) .and. n > 0))
+            end do
+          end do
+        case(2)
+          do i = 1, size(x, 2)
+            center(:, i) = merge( x(:, i) - mean_,&
+                0._dp,&
+                mask(:, i))
+          end do
+            res = matmul( center, transpose(center))
+          do j = 1, size(res, 2)
+            do i = 1, size(res, 1)
+              n = count(merge(.true., .false., mask(i, :) .and. mask(j, :)))
+              res(i, j) = res(i, j) / (n - merge(1, 0,&
+                            optval(corrected, .true.) .and. n > 0))
+            end do
+          end do
+        case default
+          call error_stop("ERROR (cov): 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 b79e39bf6..ce891bc4a 100644
--- a/src/tests/stats/CMakeLists.txt
+++ b/src/tests/stats/CMakeLists.txt
@@ -1,3 +1,4 @@
+ADDTEST(cov)
 ADDTEST(mean)
 ADDTEST(moment)
 ADDTEST(rawmoment)
diff --git a/src/tests/stats/test_cov.f90 b/src/tests/stats/test_cov.f90
new file mode 100644
index 000000000..f126e1378
--- /dev/null
+++ b/src/tests/stats/test_cov.f90
@@ -0,0 +1,566 @@
+program test_cov
+    use stdlib_experimental_error, only: check
+    use stdlib_experimental_kinds, only: sp, dp, int32, int64
+    use stdlib_experimental_stats, only: cov, var
+    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, 7._dp,&
+                                   2._dp, 4._dp, 6._dp, 8._dp,&
+                                   9._dp, 10._dp, 11._dp, 12._dp], [4, 3])
+
+    complex(dp) :: cd1(5) = [ cmplx(0.57706_dp, 0.00000_dp),&
+                            cmplx(0.00000_dp, 1.44065_dp),&
+                            cmplx(1.26401_dp, 0.00000_dp),&
+                            cmplx(0.00000_dp, 0.88833_dp),&
+                            cmplx(1.14352_dp, 0.00000_dp)]
+    complex(dp) :: ds(2,3) = reshape([ cmplx(1._dp, 0._dp),&
+                            cmplx(0._dp, 2._dp),&
+                            cmplx(3._dp, 0._dp),&
+                            cmplx(0._dp, 4._dp),&
+                            cmplx(5._dp, 0._dp),&
+                            cmplx(0._dp, 6._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(cov(x, 1) - 2.5_sp) < sptol&
+            , 'sp check 1')
+        call check( ieee_is_nan(cov(x, 1, .false.))&
+            , 'sp check 2')
+        call check( ieee_is_nan((cov(x, 1, x == 1.)))&
+            , 'sp check 3')
+        call check( abs(cov(x, 1, x < 5) - 5._sp/3) < sptol&
+            , 'sp check 4')
+        call check( abs(cov(x, 1, x < 5, corrected = .false.) -&
+            5._sp/4) < sptol&
+            , 'sp check 5')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
+            , 'sp check 6')
+        call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
+            , 'sp check 7')
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            60._sp/9, 60._sp/9, 30._sp/9&
+            ,60._sp/9, 60._sp/9, 30._sp/9&
+            ,30._sp/9, 30._sp/9, 15._sp/9]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 8')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            19._sp, 16.5_sp, 14._sp, 11.5_sp, 16.5_sp, 129._sp/9&
+            ,109.5_sp/9, 10._sp, 14._sp, 109.5_sp/9, 93._sp/9&
+            , 8.5_sp, 11.5_sp, 10._sp, 8.5_sp, 7._sp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 9')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            60._sp/9, 60._sp/9, 30._sp/9&
+            ,60._sp/9, 60._sp/9, 30._sp/9&
+            ,30._sp/9, 30._sp/9, 15._sp/9]&
+            *(size(x2, 1)-1._sp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 10')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            19._sp, 16.5_sp, 14._sp, 11.5_sp, 16.5_sp, 129._sp/9&
+            ,109.5_sp/9, 10._sp, 14._sp, 109.5_sp/9, 93._sp/9&
+            , 8.5_sp, 11.5_sp, 10._sp, 8.5_sp, 7._sp]&
+            *(size(x2, 2)-1._sp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 11')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = x2 < 10)))&
+            , 'sp check 12')
+        call check( all( abs( cov(x2, 1, mask = x2 < 11) - reshape([&
+             60._sp/9, 60._sp/9, 1._sp, 60._sp/9, 60._sp/9, 1._sp&
+             , 1._sp, 1._sp, 0.5_sp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 13')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11) - reshape([&
+            19._sp, 16.5_sp, 0.5_sp, 0.5_sp, 16.5_sp&
+            ,129._sp/9, 0.5_sp, 0.5_sp, 0.5_sp, 0.5_sp&
+            ,0.5_sp, 0.5_sp, 0.5_sp, 0.5_sp, 0.5_sp, 0.5_sp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 14')
+
+        call check( all( abs( cov(x2, 1, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+             5._sp, 5._sp, 0.5_sp, 5._sp, 5._sp, 0.5_sp, 0.5_sp,&
+             0.5_sp, 0.25_sp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'sp check 15')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+            114._sp/9, 11._sp, 0.25_sp, 0.25_sp, 11._sp, 86._sp/9,&
+            0.25_sp, 0.25_sp, 0.25_sp, 0.25_sp, 0.25_sp, 0.25_sp,&
+            0.25_sp, 0.25_sp, 0.25_sp, 0.25_sp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'sp check 16')
+
+    end subroutine test_sp
+
+    subroutine test_dp(x, x2)
+        real(dp), intent(in) :: x(:)
+        real(dp), intent(in) :: x2(:, :)
+
+        call check( abs(cov(x, 1) - 2.5_dp) < dptol&
+            , 'dp check 1')
+        call check( ieee_is_nan(cov(x, 1, .false.))&
+            , 'dp check 2')
+        call check( ieee_is_nan((cov(x, 1, x == 1.)))&
+            , 'dp check 3')
+        call check( abs(cov(x, 1, x < 5) - 5._dp/3) < dptol&
+            , 'dp check 4')
+        call check( abs(cov(x, 1, x < 5, corrected = .false.) -&
+            5._dp/4) < dptol&
+            , 'dp check 5')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
+            , 'dp check 6')
+        call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
+            , 'dp check 7')
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 8')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 9')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            *(size(x2, 1)-1._dp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 10')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            *(size(x2, 2)-1._dp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 11')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = x2 < 10)))&
+            , 'dp check 12')
+        call check( all( abs( cov(x2, 1, mask = x2 < 11) - reshape([&
+             60._dp/9, 60._dp/9, 1._dp, 60._dp/9, 60._dp/9, 1._dp&
+             , 1._dp, 1._dp, 0.5_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 13')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11) - reshape([&
+            19._dp, 16.5_dp, 0.5_dp, 0.5_dp, 16.5_dp&
+            ,129._dp/9, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp&
+            ,0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 14')
+
+        call check( all( abs( cov(x2, 1, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+             5._dp, 5._dp, 0.5_dp, 5._dp, 5._dp, 0.5_dp, 0.5_dp,&
+             0.5_dp, 0.25_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'dp check 15')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+            114._dp/9, 11._dp, 0.25_dp, 0.25_dp, 11._dp, 86._dp/9,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'dp check 16')
+
+
+    end subroutine test_dp
+
+    subroutine test_int32(x, x2)
+        integer(int32), intent(in) :: x(:)
+        integer(int32), intent(in) :: x2(:, :)
+
+        call check( abs(cov(x, 1) - 2.5_dp) < dptol&
+            , 'int32 check 1')
+        call check( ieee_is_nan(cov(x, 1, .false.))&
+            , 'int32 check 2')
+        call check( ieee_is_nan((cov(x, 1, x == 1.)))&
+            , 'int32 check 3')
+        call check( abs(cov(x, 1, x < 5) - 5._dp/3) < dptol&
+            , 'int32 check 4')
+        call check( abs(cov(x, 1, x < 5, corrected = .false.) -&
+            5._dp/4) < dptol&
+            , 'int32 check 5')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
+            , 'int32 check 6')
+        call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
+            , 'int32 check 7')
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 8')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 9')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            *(size(x2, 1)-1._dp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 10')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            *(size(x2, 2)-1._dp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 11')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = x2 < 10)))&
+            , 'int32 check 12')
+        call check( all( abs( cov(x2, 1, mask = x2 < 11) - reshape([&
+             60._dp/9, 60._dp/9, 1._dp, 60._dp/9, 60._dp/9, 1._dp&
+             , 1._dp, 1._dp, 0.5_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 13')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11) - reshape([&
+            19._dp, 16.5_dp, 0.5_dp, 0.5_dp, 16.5_dp&
+            ,129._dp/9, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp&
+            ,0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 14')
+
+        call check( all( abs( cov(x2, 1, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+             5._dp, 5._dp, 0.5_dp, 5._dp, 5._dp, 0.5_dp, 0.5_dp,&
+             0.5_dp, 0.25_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int32 check 15')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+            114._dp/9, 11._dp, 0.25_dp, 0.25_dp, 11._dp, 86._dp/9,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int32 check 16')
+
+    end subroutine test_int32
+
+    subroutine test_int64(x, x2)
+        integer(int64), intent(in) :: x(:)
+        integer(int64), intent(in) :: x2(:, :)
+
+        call check( abs(cov(x, 1) - 2.5_dp) < dptol&
+            , 'int64 check 1')
+        call check( ieee_is_nan(cov(x, 1, .false.))&
+            , 'int64 check 2')
+        call check( ieee_is_nan((cov(x, 1, x == 1)))&
+            , 'int64 check 3')
+        call check( abs(cov(x, 1, x < 5) - 5._dp/3) < dptol&
+            , 'int64 check 4')
+        call check( abs(cov(x, 1, x < 5, corrected = .false.) -&
+            5._dp/4) < dptol&
+            , 'int64 check 5')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = .false.)))&
+            , 'int64 check 6')
+        call check( any(ieee_is_nan(cov(x2, 2, mask = .false.)))&
+            , 'int64 check 7')
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 8')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 9')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            60._dp/9, 60._dp/9, 30._dp/9&
+            ,60._dp/9, 60._dp/9, 30._dp/9&
+            ,30._dp/9, 30._dp/9, 15._dp/9]&
+            *(size(x2, 1)-1._dp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 10')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            19._dp, 16.5_dp, 14._dp, 11.5_dp, 16.5_dp, 129._dp/9&
+            ,109.5_dp/9, 10._dp, 14._dp, 109.5_dp/9, 93._dp/9&
+            , 8.5_dp, 11.5_dp, 10._dp, 8.5_dp, 7._dp]&
+            *(size(x2, 2)-1._dp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 11')
+
+        call check( any(ieee_is_nan(cov(x2, 1, mask = x2 < 10)))&
+            , 'int64 check 12')
+        call check( all( abs( cov(x2, 1, mask = x2 < 11) - reshape([&
+             60._dp/9, 60._dp/9, 1._dp, 60._dp/9, 60._dp/9, 1._dp&
+             , 1._dp, 1._dp, 0.5_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 13')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11) - reshape([&
+            19._dp, 16.5_dp, 0.5_dp, 0.5_dp, 16.5_dp&
+            ,129._dp/9, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp&
+            ,0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp, 0.5_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 14')
+
+        call check( all( abs( cov(x2, 1, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+             5._dp, 5._dp, 0.5_dp, 5._dp, 5._dp, 0.5_dp, 0.5_dp,&
+             0.5_dp, 0.25_dp]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'int64 check 15')
+        call check( all( abs( cov(x2, 2, mask = x2 < 11, corrected = .false.) -&
+            reshape([&
+            114._dp/9, 11._dp, 0.25_dp, 0.25_dp, 11._dp, 86._dp/9,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp,&
+            0.25_dp, 0.25_dp, 0.25_dp, 0.25_dp]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'int64 check 16')
+
+    end subroutine test_int64
+
+    subroutine test_csp(x, x2)
+        complex(sp), intent(in) :: x(:)
+        complex(sp), intent(in) :: x2(:, :)
+
+!        complex(sp), allocatable :: cd(:,:)
+
+        call check( abs(cov(x, dim=1) -&
+                     (var(real(x),1) + var(aimag(x), 1)) ) < sptol&
+            , 'csp check 1')
+        call check( abs(cov(x, 1, aimag(x) == 0) -&
+                     var(real(x), 1, aimag(x) == 0)) < sptol&
+            , 'csp check 2')
+
+        call check( abs(cov(x, dim=1, corrected=.false.) -&
+                     (var(real(x), dim=1, corrected=.false.) +&
+                      var(aimag(x), dim=1, corrected=.false.))) <&
+                      sptol&
+            , 'csp check 3')
+
+        call check( ieee_is_nan(real(cov(x, 1, .false., corrected=.false.)))&
+            , 'csp check 4')
+
+        call check( abs(cov(x, 1, aimag(x) == 0, corrected=.false.) -&
+                      var(real(x), 1, aimag(x) == 0,&
+                      corrected=.false.)) < sptol&
+            , 'csp check 5')
+
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            (2.5_sp,0._sp), (5.5_sp,-1._sp), (8.5_sp,-2._sp)&
+            , (5.5_sp,1._sp),  (12.5_sp,0._sp), (19.5_sp,-1._sp)&
+            , (8.5_sp,2._sp),  (19.5_sp,1._sp),  (30.5_sp,0._sp)]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'csp check 6')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            (4._sp,0._sp), (0._sp,4._sp),&
+            (0._sp,-4._sp), (4._sp,0._sp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'csp check 7')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            (2.5_sp,0._sp), (5.5_sp,-1._sp), (8.5_sp,-2._sp)&
+            , (5.5_sp,1._sp),  (12.5_sp,0._sp), (19.5_sp,-1._sp)&
+            , (8.5_sp,2._sp),  (19.5_sp,1._sp),  (30.5_sp,0._sp)]&
+            *(size(x2, 1)-1._sp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < sptol)&
+            , 'csp check 8')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            (4._sp,0._sp), (0._sp,4._sp),&
+            (0._sp,-4._sp), (4._sp,0._sp)]&
+            *(size(x2, 2)-1._sp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'csp check 9')
+
+! Issue with gfortran 7 and 8: do not extract cd(1:2, 1:2) correctly
+!        allocate(cd, source = cov(x2, 1, mask = aimag(x2) < 6))
+!        call check( all( abs( cd(1:2, 1:2) - reshape([&
+!             (2.5_sp,0._sp), (5.5_sp,-1._sp)&
+!             ,(5.5_sp,1._sp), (12.5_sp,0._sp)]&
+!            ,[2, 2])&
+!            ) < sptol)&
+!            , 'csp check 10')
+!        call check( ieee_is_nan(real(cd(3,3)))&
+!            , 'csp check 10 bis')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6) - reshape([&
+             (4._sp,0._sp), (0._sp,2._sp)&
+             ,(0._sp,-2._sp), (2._sp,0._sp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sptol)&
+            , 'csp check 11')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6, corrected = .false.) -&
+            reshape([&
+             (2.6666666666666666_sp,0._sp), (0._sp,1._sp)&
+             ,(0._sp,-1._sp), (1._sp,0._sp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < sp)&
+            , 'csp check 12')
+
+    end subroutine test_csp
+
+    subroutine test_cdp(x, x2)
+        complex(dp), intent(in) :: x(:)
+        complex(dp), intent(in) :: x2(:, :)
+
+!        complex(dp), allocatable :: cd(:,:)
+
+        call check( abs(cov(x, dim=1) -&
+                     (var(real(x),1) + var(aimag(x), 1)) ) < dptol&
+            , 'cdp check 1')
+        call check( abs(cov(x, 1, aimag(x) == 0) -&
+                     var(real(x), 1, aimag(x) == 0)) < dptol&
+            , 'cdp check 2')
+
+        call check( abs(cov(x, dim=1, corrected=.false.) -&
+                     (var(real(x), dim=1, corrected=.false.) +&
+                      var(aimag(x), dim=1, corrected=.false.))) <&
+                      dptol&
+            , 'cdp check 3')
+
+        call check( ieee_is_nan(real(cov(x, 1, .false., corrected=.false.)))&
+            , 'cdp check 4')
+
+        call check( abs(cov(x, 1, aimag(x) == 0, corrected=.false.) -&
+                      var(real(x), 1, aimag(x) == 0,&
+                      corrected=.false.)) < dptol&
+            , 'cdp check 5')
+
+
+        call check( all( abs( cov(x2, 1) - reshape([&
+            (2.5_dp,0._dp), (5.5_dp,-1._dp), (8.5_dp,-2._dp)&
+            , (5.5_dp,1._dp),  (12.5_dp,0._dp), (19.5_dp,-1._dp)&
+            , (8.5_dp,2._dp),  (19.5_dp,1._dp),  (30.5_dp,0._dp)]&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'cdp check 6')
+        call check( all( abs( cov(x2, 2) - reshape([&
+            (4._dp,0._dp), (0._dp,4._dp),&
+            (0._dp,-4._dp), (4._dp,0._dp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 7')
+
+        call check( all( abs( cov(x2, 1, corrected=.false.) - reshape([&
+            (2.5_dp,0._dp), (5.5_dp,-1._dp), (8.5_dp,-2._dp)&
+            , (5.5_dp,1._dp),  (12.5_dp,0._dp), (19.5_dp,-1._dp)&
+            , (8.5_dp,2._dp),  (19.5_dp,1._dp),  (30.5_dp,0._dp)]&
+            *(size(x2, 1)-1._dp)/size(x2, 1)&
+            ,[ size(x2, 2), size(x2, 2)])&
+            ) < dptol)&
+            , 'cdp check 8')
+        call check( all( abs( cov(x2, 2, corrected=.false.) - reshape([&
+            (4._dp,0._dp), (0._dp,4._dp),&
+            (0._dp,-4._dp), (4._dp,0._dp)]&
+            *(size(x2, 2)-1._dp)/size(x2, 2)&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 9')
+
+! Issue with gfortran 7 and 8: do not extract cd(1:2, 1:2) correctly
+!        allocate(cd, source = cov(x2, 1, mask = aimag(x2) < 6))
+!
+!        call check( all( abs( cd(1:2, 1:2) - reshape([&
+!             (2.5_dp,0._dp), (5.5_dp,-1._dp)&
+!             ,(5.5_dp,1._dp), (12.5_dp,0._dp)]&
+!            ,[2, 2])&
+!            ) < dptol)&
+!            , 'cdp check 10')
+!        call check( ieee_is_nan(real(cd(3,3)))&
+!            , 'cdp check 10 bis')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6) - reshape([&
+             (4._dp,0._dp), (0._dp,2._dp)&
+             ,(0._dp,-2._dp), (2._dp,0._dp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 11')
+
+        call check( all( abs( cov(x2, 2, mask = aimag(x2) < 6, corrected = .false.) -&
+            reshape([&
+             (2.6666666666666666_dp,0._dp), (0._dp,1._dp)&
+             ,(0._dp,-1._dp), (1._dp,0._dp)]&
+            ,[ size(x2, 1), size(x2, 1)])&
+            ) < dptol)&
+            , 'cdp check 12')
+
+    end subroutine test_cdp
+
+end program test_cov

From 7b1ebc799eab44bc962afe4bf4b779d4b71b1ee0 Mon Sep 17 00:00:00 2001
From: "Vandenplas, Jeremie" <jeremie.vandenplas@gmail.com>
Date: Sun, 12 Apr 2020 11:19:12 +0200
Subject: [PATCH 2/3] covariance_dev: spec

---
 src/stdlib_experimental_stats.md | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md
index 671351715..c4597b5da 100644
--- a/src/stdlib_experimental_stats.md
+++ b/src/stdlib_experimental_stats.md
@@ -10,12 +10,12 @@
 
 ### Description
 
-Returns the covariance of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true` (if `mask` is provided).
+Returns the covariance of the elements of `array` along dimension `dim` if the corresponding element in `mask` is `true`.
 
 Per default, the covariance is defined as:
 
 ```
- cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)) )
+ cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)))
 ```
 
 where n is the number of elements.
@@ -29,7 +29,7 @@ The scaling can be changed with the logical argument `corrected`. If `corrected`
 
 ### Arguments
 
-`array`: Shall be a 1-rank or a 2-rank array of type `integer`, `real`, or `complex`.
+`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`.
 
@@ -39,14 +39,15 @@ The scaling can be changed with the logical argument `corrected`. If `corrected`
 
 ### Return value
 
-If `array` is of rank 1 and of type `real` or `complex`, the result is of the same type `real` corresponding to the type of `array`.
+If `array` is of rank 1 and of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
 If `array` is of rank 2 and of type `real` or `complex`, the result is of the same type as `array`.
 If `array` is of type `integer`, the result is of type `real(dp)`.
 
 If `array` is of rank 1, a scalar with the covariance (that is the variance) of all elements in `array` is returned.
-If `array` is of rank 2, a 2-rank array is returned.
+If `array` is of rank 2, a rank-2 array is returned.
 
 If `mask` is specified, the result is the covariance 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

From 8f626cdcd73525c03ea31a9abb059c2329251e1c Mon Sep 17 00:00:00 2001
From: "Vandenplas, Jeremie" <jeremie.vandenplas@gmail.com>
Date: Mon, 13 Apr 2020 19:46:06 +0200
Subject: [PATCH 3/3] covariance_dev: following @ivan-pi s comments

---
 src/stdlib_experimental_stats.md | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/stdlib_experimental_stats.md b/src/stdlib_experimental_stats.md
index c4597b5da..8fe1d763f 100644
--- a/src/stdlib_experimental_stats.md
+++ b/src/stdlib_experimental_stats.md
@@ -18,7 +18,7 @@ Per default, the covariance is defined as:
  cov(array) = 1/(n-1) sum_i (array(i) - mean(array) * (array(i) - mean(array)))
 ```
 
-where n is the number of elements.
+where `n` is the number of elements.
 
 The scaling can be changed with the logical argument `corrected`. If `corrected` is `.false.`, then the sum is scaled with `n`, otherwise with `n-1`.
 
@@ -31,9 +31,9 @@ The scaling can be changed with the logical argument `corrected`. If `corrected`
 
 `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`.
+`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 by a scalar or an array of the same shape as `array`.
+`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
 
 `corrected` (optional): Shall be a scalar of type `logical`. If `corrected` is `.true.` (default value), the sum is scaled with `n-1`. If `corrected` is `.false.`, then the sum is scaled with `n`.
 
@@ -78,16 +78,16 @@ Returns the mean of all the elements of `array`, or of the elements of `array` a
 
 `array`: Shall be an 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`.
+`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 by a scalar or an array of the same shape as `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 type `real` or `complex`, the result is of the same type as `array`.
 If `array` is of type `integer`, the result is of type `real(dp)`.
 
-If `dim` is absent, a scalar with the mean of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
+If `dim` is absent, a scalar with the mean of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
 
 If `mask` is specified, the result is the mean of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
 
@@ -121,7 +121,7 @@ The _k_-th order central moment is defined as :
  moment(array) = 1/n sum_i (array(i) - mean(array))^k
 ```
 
-where n is the number of elements.
+where `n` is the number of elements.
 
 The _k_-th order moment about `center` is defined as :
 
@@ -141,18 +141,18 @@ The _k_-th order moment about `center` is defined as :
 
 `order`: Shall be an scalar of type `integer`.
 
-`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`.
+`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`.
 
 `center` (optional): Shall be a scalar of the same type of `result` if `dim` is not provided. If `dim` is provided, `center` shall be a scalar or an array (with a shape similar to that of `array` with dimension `dim` dropped) of the same type of `result`.
 
-`mask` (optional): Shall be of type `logical` and either by a scalar or an array of the same shape as `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 type `real` or `complex`, the result is of the same type as `array`.
 If `array` is of type `integer`, the result is of type `real(dp)`.
 
-If `dim` is absent, a scalar with the _k_-th (central) moment of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
+If `dim` is absent, a scalar with the _k_-th (central) moment of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
 
 If `mask` is specified, the result is the _k_-th  (central) moment of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.
 
@@ -185,7 +185,7 @@ Per default, the variance is defined as the best unbiased estimator and is compu
  var(array) = 1/(n-1) sum_i (array(i) - mean(array))^2
 ```
 
-where n is the number of elements.
+where `n` is the number of elements.
 
 The use of the term `n-1` for scaling is called Bessel 's correction. The scaling can be changed with the logical argument `corrected`. If `corrected` is `.false.`, then the sum is scaled with `n`, otherwise with `n-1`.
 
@@ -200,9 +200,9 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
 
 `array`: Shall be an 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`.
+`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 by a scalar or an array of the same shape as `array`.
+`mask` (optional): Shall be of type `logical` and either a scalar or an array of the same shape as `array`.
 
 `corrected` (optional): Shall be a scalar of type `logical`. If `corrected` is `.true.` (default value), the sum is scaled with `n-1`. If `corrected` is `.false.`, then the sum is scaled with `n`.
 
@@ -211,7 +211,7 @@ The use of the term `n-1` for scaling is called Bessel 's correction. The scalin
 If `array` is of type `real` or `complex`, the result is of type `real` corresponding to the type of `array`.
 If `array` is of type `integer`, the result is of type `real(dp)`.
 
-If `dim` is absent, a scalar with the variance of all elements in `array` is returned. Otherwise, an array of rank n-1, where n equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
+If `dim` is absent, a scalar with the variance of all elements in `array` is returned. Otherwise, an array of rank `n-1`, where `n` equals the rank of `array`, and a shape similar to that of `array` with dimension `dim` dropped is returned.
 
 If `mask` is specified, the result is the variance of all elements of `array` corresponding to `true` elements of `mask`. If every element of `mask` is `false`, the result is IEEE `NaN`.