From 02a0f011d4b4fb798df62afb0244c0d7d9f7a107 Mon Sep 17 00:00:00 2001
From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
Date: Mon, 28 Dec 2020 10:59:30 +0100
Subject: [PATCH 1/2] Reduce compilation load by splitting stdlib_stats_moment

- max rank 15 and one ninja job result in ~7 GB peak memory usage
---
 src/CMakeLists.txt                  |   2 +
 src/Makefile.manual                 |   4 +
 src/stdlib_stats_moment.fypp        | 208 ----------------------------
 src/stdlib_stats_moment_all.fypp    | 123 ++++++++++++++++
 src/stdlib_stats_moment_scalar.fypp | 110 +++++++++++++++
 5 files changed, 239 insertions(+), 208 deletions(-)
 create mode 100644 src/stdlib_stats_moment_all.fypp
 create mode 100644 src/stdlib_stats_moment_scalar.fypp

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 02604959e..a57e2731b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -14,6 +14,8 @@ set(fppFiles
     stdlib_stats_cov.fypp
     stdlib_stats_mean.fypp
     stdlib_stats_moment.fypp
+    stdlib_stats_moment_all.fypp
+    stdlib_stats_moment_scalar.fypp
     stdlib_stats_var.fypp
     stdlib_quadrature.fypp
     stdlib_quadrature_trapz.fypp
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 872f704c0..c300570ef 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -15,6 +15,8 @@ SRC = f18estop.f90 \
       stdlib_stats.f90 \
       stdlib_stats_mean.f90 \
       stdlib_stats_moment.f90 \
+      stdlib_stats_moment_all.f90 \
+      stdlib_stats_moment_scalar.f90 \
       stdlib_stats_var.f90
 
 LIB = libstdlib.a
@@ -79,4 +81,6 @@ stdlib_quadrature.f90: stdlib_quadrature.fypp
 stdlib_stats.f90: stdlib_stats.fypp
 stdlib_stats_mean.f90: stdlib_stats_mean.fypp
 stdlib_stats_moment.f90: stdlib_stats_moment.fypp
+stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp
+stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp
 stdlib_stats_var.f90: stdlib_stats_var.fypp
diff --git a/src/stdlib_stats_moment.fypp b/src/stdlib_stats_moment.fypp
index ab1b77a5d..578a86aea 100644
--- a/src/stdlib_stats_moment.fypp
+++ b/src/stdlib_stats_moment.fypp
@@ -11,93 +11,6 @@ submodule (stdlib_stats) stdlib_stats_moment
 
 contains
 
-  #:for k1, t1 in RC_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_all",rank, t1, k1)
-      module function ${RName}$(x, order, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        ${t1}$, intent(in), optional :: center
-        logical, intent(in), optional :: mask
-        ${t1}$ :: res
-
-        real(${k1}$) :: n
-
-        if (.not.optval(mask, .true.)) then
-          res = ieee_value(1._${k1}$, ieee_quiet_nan)
-          return
-        end if
-
-        n = real(size(x, kind = int64), ${k1}$)
-
-        if (present(center)) then
-         res = sum((x - center)**order) / n
-        else
-         res = sum((x - mean(x))**order) / n
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
-  #:for k1, t1 in INT_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_all",rank, t1, k1, 'dp')
-      module function ${RName}$(x, order, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        real(dp), intent(in), optional :: center
-        logical, intent(in), optional :: mask
-        real(dp) :: res
-
-        real(dp) :: n
-
-        if (.not.optval(mask, .true.)) then
-          res = ieee_value(1._dp, ieee_quiet_nan)
-          return
-        end if
-
-        n = real(size(x, kind = int64), dp)
-
-        if (present(center)) then
-         res = sum((real(x, dp) - center)**order) / n
-        else
-         res = sum((real(x, dp) - mean(x))**order) / n
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
-  #:for k1, t1 in RC_KINDS_TYPES
-    #:for rank in REDRANKS
-      #:set RName = rname("moment_scalar",rank, t1, k1)
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        ${t1}$, intent(in) :: center
-        logical, intent(in), optional :: mask
-        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
-
-        if (.not.optval(mask, .true.)) then
-          res = ieee_value(1._${k1}$, ieee_quiet_nan)
-          return
-        end if
-
-        if (dim >= 1 .and. dim <= ${rank}$) then
-          res = sum((x - center)**order, dim) / size(x, dim)
-        else
-          call error_stop("ERROR (moment): wrong dimension")
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("moment",rank, t1, k1)
@@ -146,33 +59,6 @@ contains
   #:endfor
 
 
-  #:for k1, t1 in INT_KINDS_TYPES
-    #:for rank in REDRANKS
-      #:set RName = rname("moment_scalar",rank, t1, k1, 'dp')
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        real(dp),intent(in) :: center
-        logical, intent(in), optional :: mask
-        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
-
-        if (.not.optval(mask, .true.)) then
-          res = ieee_value(1._dp, ieee_quiet_nan)
-          return
-        end if
-
-        if (dim >= 1 .and. dim <= ${rank}$) then
-          res = sum( (real(x, dp) - center)**order, dim) / size(x, dim)
-        else
-          call error_stop("ERROR (moment): wrong dimension")
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("moment",rank, t1, k1, 'dp')
@@ -222,78 +108,6 @@ contains
   #:endfor
 
 
-  #:for k1, t1 in RC_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_mask_all",rank, t1, k1)
-      module function ${RName}$(x, order, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        ${t1}$, intent(in), optional :: center
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        ${t1}$ :: res
-
-        real(${k1}$) :: n
-
-        n = real(count(mask, kind = int64), ${k1}$)
-
-        if (present(center)) then
-         res = sum((x - center)**order, mask) / n
-        else
-         res = sum((x - mean(x, mask))**order, mask) / n
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
-  #:for k1, t1 in INT_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp')
-      module function ${RName}$(x, order, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        real(dp),intent(in), optional :: center
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        real(dp) :: res
-
-        real(dp) :: n
-
-        n = real(count(mask, kind = int64), dp)
-
-        if (present(center)) then
-         res = sum((real(x, dp) - center)**order, mask) / n
-        else
-         res = sum((real(x, dp) - mean(x,mask))**order, mask) / n
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
-  #:for k1, t1 in RC_KINDS_TYPES
-    #:for rank in REDRANKS
-      #:set RName = rname("moment_mask_scalar",rank, t1, k1)
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        ${t1}$, intent(in) :: center
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
-
-        if (dim >= 1 .and. dim <= ${rank}$) then
-          res = sum((x - center)**order, dim, mask) / count(mask, dim)
-        else
-          call error_stop("ERROR (moment): wrong dimension")
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
   #:for k1, t1 in RC_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("moment_mask",rank, t1, k1)
@@ -350,28 +164,6 @@ contains
   #:endfor
 
 
-  #:for k1, t1 in INT_KINDS_TYPES
-    #:for rank in REDRANKS
-      #:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp')
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        real(dp), intent(in) :: center
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
-
-        if (dim >= 1 .and. dim <= ${rank}$) then
-          res = sum(( real(x, dp) - center)**order, dim, mask) / count(mask, dim)
-        else
-          call error_stop("ERROR (moment): wrong dimension")
-        end if
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
   #:for k1, t1 in INT_KINDS_TYPES
     #:for rank in RANKS
       #:set RName = rname("moment_mask",rank, t1, k1, 'dp')
diff --git a/src/stdlib_stats_moment_all.fypp b/src/stdlib_stats_moment_all.fypp
new file mode 100644
index 000000000..34aef0d56
--- /dev/null
+++ b/src/stdlib_stats_moment_all.fypp
@@ -0,0 +1,123 @@
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set REDRANKS = range(2, MAXRANK + 1)
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_stats) stdlib_stats_moment_all
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+  use stdlib_error, only: error_stop
+  use stdlib_optval, only: optval
+  implicit none
+
+contains
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_all",rank, t1, k1)
+      module function ${RName}$(x, order, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        ${t1}$, intent(in), optional :: center
+        logical, intent(in), optional :: mask
+        ${t1}$ :: res
+
+        real(${k1}$) :: n
+
+        if (.not.optval(mask, .true.)) then
+          res = ieee_value(1._${k1}$, ieee_quiet_nan)
+          return
+        end if
+
+        n = real(size(x, kind = int64), ${k1}$)
+
+        if (present(center)) then
+         res = sum((x - center)**order) / n
+        else
+         res = sum((x - mean(x))**order) / n
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_all",rank, t1, k1, 'dp')
+      module function ${RName}$(x, order, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        real(dp), intent(in), optional :: center
+        logical, intent(in), optional :: mask
+        real(dp) :: res
+
+        real(dp) :: n
+
+        if (.not.optval(mask, .true.)) then
+          res = ieee_value(1._dp, ieee_quiet_nan)
+          return
+        end if
+
+        n = real(size(x, kind = int64), dp)
+
+        if (present(center)) then
+         res = sum((real(x, dp) - center)**order) / n
+        else
+         res = sum((real(x, dp) - mean(x))**order) / n
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_mask_all",rank, t1, k1)
+      module function ${RName}$(x, order, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        ${t1}$, intent(in), optional :: center
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        ${t1}$ :: res
+
+        real(${k1}$) :: n
+
+        n = real(count(mask, kind = int64), ${k1}$)
+
+        if (present(center)) then
+         res = sum((x - center)**order, mask) / n
+        else
+         res = sum((x - mean(x, mask))**order, mask) / n
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_mask_all",rank, t1, k1, 'dp')
+      module function ${RName}$(x, order, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        real(dp),intent(in), optional :: center
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        real(dp) :: res
+
+        real(dp) :: n
+
+        n = real(count(mask, kind = int64), dp)
+
+        if (present(center)) then
+         res = sum((real(x, dp) - center)**order, mask) / n
+        else
+         res = sum((real(x, dp) - mean(x,mask))**order, mask) / n
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+end submodule
diff --git a/src/stdlib_stats_moment_scalar.fypp b/src/stdlib_stats_moment_scalar.fypp
new file mode 100644
index 000000000..137ff6d79
--- /dev/null
+++ b/src/stdlib_stats_moment_scalar.fypp
@@ -0,0 +1,110 @@
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set REDRANKS = range(2, MAXRANK + 1)
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_stats) stdlib_stats_moment_scalar
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+  use stdlib_error, only: error_stop
+  use stdlib_optval, only: optval
+  implicit none
+
+contains
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:for rank in REDRANKS
+      #:set RName = rname("moment_scalar",rank, t1, k1)
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        ${t1}$, intent(in) :: center
+        logical, intent(in), optional :: mask
+        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
+
+        if (.not.optval(mask, .true.)) then
+          res = ieee_value(1._${k1}$, ieee_quiet_nan)
+          return
+        end if
+
+        if (dim >= 1 .and. dim <= ${rank}$) then
+          res = sum((x - center)**order, dim) / size(x, dim)
+        else
+          call error_stop("ERROR (moment): wrong dimension")
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:for rank in REDRANKS
+      #:set RName = rname("moment_scalar",rank, t1, k1, 'dp')
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        real(dp),intent(in) :: center
+        logical, intent(in), optional :: mask
+        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
+
+        if (.not.optval(mask, .true.)) then
+          res = ieee_value(1._dp, ieee_quiet_nan)
+          return
+        end if
+
+        if (dim >= 1 .and. dim <= ${rank}$) then
+          res = sum( (real(x, dp) - center)**order, dim) / size(x, dim)
+        else
+          call error_stop("ERROR (moment): wrong dimension")
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:for rank in REDRANKS
+      #:set RName = rname("moment_mask_scalar",rank, t1, k1)
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        ${t1}$, intent(in) :: center
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
+
+        if (dim >= 1 .and. dim <= ${rank}$) then
+          res = sum((x - center)**order, dim, mask) / count(mask, dim)
+        else
+          call error_stop("ERROR (moment): wrong dimension")
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:for rank in REDRANKS
+      #:set RName = rname("moment_mask_scalar",rank, t1, k1, 'dp')
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        real(dp), intent(in) :: center
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
+
+        if (dim >= 1 .and. dim <= ${rank}$) then
+          res = sum(( real(x, dp) - center)**order, dim, mask) / count(mask, dim)
+        else
+          call error_stop("ERROR (moment): wrong dimension")
+        end if
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+end submodule

From 27fe4219552067137be02d38eca3d3afc5e3d743 Mon Sep 17 00:00:00 2001
From: Sebastian Ehlert <28669218+awvwgk@users.noreply.github.com>
Date: Mon, 28 Dec 2020 11:08:23 +0100
Subject: [PATCH 2/2] Split stdlib_stats_moment again to separate the masked
 procedures

- with default max rank and two ninja jobs ~7 GB peak memory usage
---
 src/CMakeLists.txt                |   1 +
 src/Makefile.manual               |   2 +
 src/stdlib_stats_moment.fypp      | 103 --------------------------
 src/stdlib_stats_moment_mask.fypp | 116 ++++++++++++++++++++++++++++++
 4 files changed, 119 insertions(+), 103 deletions(-)
 create mode 100644 src/stdlib_stats_moment_mask.fypp

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index a57e2731b..e2caa0bbc 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -15,6 +15,7 @@ set(fppFiles
     stdlib_stats_mean.fypp
     stdlib_stats_moment.fypp
     stdlib_stats_moment_all.fypp
+    stdlib_stats_moment_mask.fypp
     stdlib_stats_moment_scalar.fypp
     stdlib_stats_var.fypp
     stdlib_quadrature.fypp
diff --git a/src/Makefile.manual b/src/Makefile.manual
index c300570ef..dd6f12708 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -16,6 +16,7 @@ SRC = f18estop.f90 \
       stdlib_stats_mean.f90 \
       stdlib_stats_moment.f90 \
       stdlib_stats_moment_all.f90 \
+      stdlib_stats_moment_mask.f90 \
       stdlib_stats_moment_scalar.f90 \
       stdlib_stats_var.f90
 
@@ -82,5 +83,6 @@ stdlib_stats.f90: stdlib_stats.fypp
 stdlib_stats_mean.f90: stdlib_stats_mean.fypp
 stdlib_stats_moment.f90: stdlib_stats_moment.fypp
 stdlib_stats_moment_all.f90: stdlib_stats_moment_all.fypp
+stdlib_stats_moment_mask.f90: stdlib_stats_moment_mask.fypp
 stdlib_stats_moment_scalar.f90: stdlib_stats_moment_scalar.fypp
 stdlib_stats_var.f90: stdlib_stats_var.fypp
diff --git a/src/stdlib_stats_moment.fypp b/src/stdlib_stats_moment.fypp
index 578a86aea..7283cc576 100644
--- a/src/stdlib_stats_moment.fypp
+++ b/src/stdlib_stats_moment.fypp
@@ -107,107 +107,4 @@ contains
     #:endfor
   #:endfor
 
-
-  #:for k1, t1 in RC_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_mask",rank, t1, k1)
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
-
-        integer :: i
-        real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
-        ${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$
-
-        n = real(count(mask, dim), ${k1}$)
-
-        res = 0
-        select case(dim)
-          #:for fi in range(1, rank+1)
-          case(${fi}$)
-            if (present(center)) then
-              do i = 1, size(x, ${fi}$)
-                res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ -&
-                  center)**order,&
-                  #:if t1[0] == 'r'
-                    0._${k1}$,&
-                  #:else
-                    cmplx(0,0,kind=${k1}$),&
-                  #:endif
-                    mask${select_subarray(rank, [(fi, 'i')])}$)
-              end do
-            else
-              allocate(mean_, source = mean(x, ${fi}$, mask))
-              do i = 1, size(x, ${fi}$)
-                res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order,&
-                  #:if t1[0] == 'r'
-                    0._${k1}$,&
-                  #:else
-                    cmplx(0,0,kind=${k1}$),&
-                  #:endif
-                    mask${select_subarray(rank, [(fi, 'i')])}$)
-              end do
-              deallocate(mean_)
-            end if
-          #:endfor
-          case default
-            call error_stop("ERROR (moment): wrong dimension")
-        end select
-        res = res / n
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
-
-  #:for k1, t1 in INT_KINDS_TYPES
-    #:for rank in RANKS
-      #:set RName = rname("moment_mask",rank, t1, k1, 'dp')
-      module function ${RName}$(x, order, dim, center, mask) result(res)
-        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
-        integer, intent(in) :: order
-        integer, intent(in) :: dim
-        real(dp), intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
-        logical, intent(in) :: mask${ranksuffix(rank)}$
-        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
-
-        integer :: i
-        real(dp) :: n${reduced_shape('x', rank, 'dim')}$
-        real(dp), allocatable :: mean_${ranksuffix(rank-1)}$
-
-        n = real(count(mask, dim), dp)
-
-        res = 0
-        select case(dim)
-          #:for fi in range(1, rank+1)
-          case(${fi}$)
-            if (present(center)) then
-              do i = 1, size(x, ${fi}$)
-                res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -&
-                                    center)**order,&
-                                    0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
-              end do
-            else
-              allocate(mean_, source = mean(x, ${fi}$, mask))
-              do i = 1, size(x, ${fi}$)
-                res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)&
-                                    **order,&
-                                    0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
-              end do
-              deallocate(mean_)
-            end if
-          #:endfor
-          case default
-            call error_stop("ERROR (moment): wrong dimension")
-        end select
-        res = res / n
-
-      end function ${RName}$
-    #:endfor
-  #:endfor
-
 end submodule
diff --git a/src/stdlib_stats_moment_mask.fypp b/src/stdlib_stats_moment_mask.fypp
new file mode 100644
index 000000000..a56d4d1b3
--- /dev/null
+++ b/src/stdlib_stats_moment_mask.fypp
@@ -0,0 +1,116 @@
+#:include "common.fypp"
+#:set RANKS = range(1, MAXRANK + 1)
+#:set REDRANKS = range(2, MAXRANK + 1)
+#:set RC_KINDS_TYPES = REAL_KINDS_TYPES + CMPLX_KINDS_TYPES
+submodule (stdlib_stats) stdlib_stats_moment_mask
+
+  use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
+  use stdlib_error, only: error_stop
+  use stdlib_optval, only: optval
+  implicit none
+
+contains
+
+  #:for k1, t1 in RC_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_mask",rank, t1, k1)
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        ${t1}$, intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        ${t1}$ :: res${reduced_shape('x', rank, 'dim')}$
+
+        integer :: i
+        real(${k1}$) :: n${reduced_shape('x', rank, 'dim')}$
+        ${t1}$, allocatable :: mean_${ranksuffix(rank-1)}$
+
+        n = real(count(mask, dim), ${k1}$)
+
+        res = 0
+        select case(dim)
+          #:for fi in range(1, rank+1)
+          case(${fi}$)
+            if (present(center)) then
+              do i = 1, size(x, ${fi}$)
+                res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ -&
+                  center)**order,&
+                  #:if t1[0] == 'r'
+                    0._${k1}$,&
+                  #:else
+                    cmplx(0,0,kind=${k1}$),&
+                  #:endif
+                    mask${select_subarray(rank, [(fi, 'i')])}$)
+              end do
+            else
+              allocate(mean_, source = mean(x, ${fi}$, mask))
+              do i = 1, size(x, ${fi}$)
+                res = res + merge( (x${select_subarray(rank, [(fi, 'i')])}$ - mean_)**order,&
+                  #:if t1[0] == 'r'
+                    0._${k1}$,&
+                  #:else
+                    cmplx(0,0,kind=${k1}$),&
+                  #:endif
+                    mask${select_subarray(rank, [(fi, 'i')])}$)
+              end do
+              deallocate(mean_)
+            end if
+          #:endfor
+          case default
+            call error_stop("ERROR (moment): wrong dimension")
+        end select
+        res = res / n
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+
+  #:for k1, t1 in INT_KINDS_TYPES
+    #:for rank in RANKS
+      #:set RName = rname("moment_mask",rank, t1, k1, 'dp')
+      module function ${RName}$(x, order, dim, center, mask) result(res)
+        ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+        integer, intent(in) :: order
+        integer, intent(in) :: dim
+        real(dp), intent(in), optional :: center${reduced_shape('x', rank, 'dim')}$
+        logical, intent(in) :: mask${ranksuffix(rank)}$
+        real(dp) :: res${reduced_shape('x', rank, 'dim')}$
+
+        integer :: i
+        real(dp) :: n${reduced_shape('x', rank, 'dim')}$
+        real(dp), allocatable :: mean_${ranksuffix(rank-1)}$
+
+        n = real(count(mask, dim), dp)
+
+        res = 0
+        select case(dim)
+          #:for fi in range(1, rank+1)
+          case(${fi}$)
+            if (present(center)) then
+              do i = 1, size(x, ${fi}$)
+                res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) -&
+                                    center)**order,&
+                                    0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
+              end do
+            else
+              allocate(mean_, source = mean(x, ${fi}$, mask))
+              do i = 1, size(x, ${fi}$)
+                res = res + merge((real(x${select_subarray(rank, [(fi, 'i')])}$, dp) - mean_)&
+                                    **order,&
+                                    0._dp, mask${select_subarray(rank, [(fi, 'i')])}$)
+              end do
+              deallocate(mean_)
+            end if
+          #:endfor
+          case default
+            call error_stop("ERROR (moment): wrong dimension")
+        end select
+        res = res / n
+
+      end function ${RName}$
+    #:endfor
+  #:endfor
+
+end submodule