jvdp1 · Jan 21, 2020
diff --git a/‎.github/workflows/CI.yml
+8-2 b/‎.github/workflows/CI.yml
+8-2
diff --git a/‎.github/workflows/ci_windows.yml
+3-1 b/‎.github/workflows/ci_windows.yml
+3-1
diff --git a/‎CMakeLists.txt
+1-5 b/‎CMakeLists.txt
+1-5
diff --git a/‎README.md
+26-7 b/‎README.md
+26-7
diff --git a/‎WORKFLOW.md
+55-3 b/‎WORKFLOW.md
+55-3
diff --git a/‎src/CMakeLists.txt
+19-1 b/‎src/CMakeLists.txt
+19-1
diff --git a/‎src/Makefile.manual
+7-1 b/‎src/Makefile.manual
+7-1
diff --git a/‎src/f03_stdlib_experimental_stats.f90
+1,667 b/‎src/f03_stdlib_experimental_stats.f90
+1,667
diff --git a/‎src/f03_stdlib_experimental_stats_mean.f90
+25,896 b/‎src/f03_stdlib_experimental_stats_mean.f90
+25,896
diff --git a/‎src/f90_stdlib_experimental_stats.f90
+575 b/‎src/f90_stdlib_experimental_stats.f90
+575
diff --git a/‎src/f90_stdlib_experimental_stats_mean.f90
+4,225 b/‎src/f90_stdlib_experimental_stats_mean.f90
+4,225
diff --git a/‎src/stdlib_experimental_io.f90
+7-4 b/‎src/stdlib_experimental_io.f90
+7-4
diff --git a/‎src/stdlib_experimental_optval.f90
+8-8 b/‎src/stdlib_experimental_optval.f90
+8-8
diff --git a/‎src/stdlib_experimental_stats.fypp
+124 b/‎src/stdlib_experimental_stats.fypp
+124
diff --git a/‎src/stdlib_experimental_stats.md
+43 b/‎src/stdlib_experimental_stats.md
+43
diff --git a/‎src/stdlib_experimental_stats_mean.fypp
+293 b/‎src/stdlib_experimental_stats_mean.fypp
+293
diff --git a/‎src/stdlib_experimental_system.F90
+43 b/‎src/stdlib_experimental_system.F90
+43
diff --git a/‎src/tests/CMakeLists.txt
+2 b/‎src/tests/CMakeLists.txt
+2
diff --git a/‎src/tests/Makefile.manual
+3 b/‎src/tests/Makefile.manual
+3
diff --git a/‎src/tests/optval/test_optval.f90
+139-18 b/‎src/tests/optval/test_optval.f90
+139-18
diff --git a/‎src/tests/stats/CMakeLists.txt
+5 b/‎src/tests/stats/CMakeLists.txt
+5
diff --git a/‎src/tests/stats/Makefile.manual
+3 b/‎src/tests/stats/Makefile.manual
+3
diff --git a/‎src/tests/stats/array3.dat
+16 b/‎src/tests/stats/array3.dat
+16
diff --git a/‎src/tests/stats/test_mean.f90
+75 b/‎src/tests/stats/test_mean.f90
+75
diff --git a/‎src/tests/stats/test_mean_f03.f90
+38 b/‎src/tests/stats/test_mean_f03.f90
+38
diff --git a/‎src/tests/system/CMakeLists.txt
+5 b/‎src/tests/system/CMakeLists.txt
+5
diff --git a/‎src/tests/system/test_sleep.f90
+33 b/‎src/tests/system/test_sleep.f90
+33
@@ -7,6 +7,7 @@ env:
   CMAKE_BUILD_PARALLEL_LEVEL: "2" # 2 cores on each GHA VM, enable parallel builds
   CTEST_OUTPUT_ON_FAILURE: "ON" # This way we don't need a flag to ctest
   CTEST_PARALLEL_LEVEL: "2"
+  CTEST_TIME_TIMEOUT: "5"  # some failures hang forever
   HOMEBREW_NO_ANALYTICS: "ON" # Make Homebrew installation a little quicker
   HOMEBREW_NO_AUTO_UPDATE: "ON"
   HOMEBREW_NO_BOTTLE_SOURCE_FALLBACK: "ON"
@@ -55,10 +56,15 @@ jobs:
       run: cmake -Wdev -DCMAKE_BUILD_TYPE=Release -S . -B build
 
     - name: Build and compile
-      run: cmake --build build || cmake --build build --verbose --parallel 1
+      run: cmake --build build
+
+    - name: catch build fail
+      run: cmake --build build --verbose --parallel 1
+      if: failure()
 
     - name: test
-      run: cmake --build build --target test
+      run: ctest --parallel --output-on-failure
+      working-directory: build
 
     - name: Test in-tree builds
       if: contains( matrix.gcc_v, '9') # Only test one compiler on each platform
 
@@ -4,6 +4,7 @@ on: [push, pull_request]
 
 env:
   CI: "ON"
+  CTEST_TIME_TIMEOUT: "5"  # some failures hang forever
 
 jobs:
   Build:
@@ -24,7 +25,8 @@ jobs:
     - name: CMake build
       run: cmake --build build --parallel
 
-    - run: cmake --build build --verbose --parallel 1
+    - name: catch build fail
+      run: cmake --build build --verbose --parallel 1
       if: failure()
 
     - name: CTest
 
@@ -2,11 +2,6 @@ cmake_minimum_required(VERSION 3.14.0)
 project(stdlib Fortran)
 enable_testing()
 
-# this avoids stdlib and projects using stdlib from having to introspect stdlib's directory structure
-# FIXME: this eventually needs to be handled more precisely, as this spills all .mod/.smod into one directory
-#  and thereby can clash if module/submodule names are the same in different parts of library
-set(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_BINARY_DIR})
-
 # --- compiler options
 if(CMAKE_Fortran_COMPILER_ID STREQUAL GNU)
   add_compile_options(-fimplicit-none)
@@ -20,6 +15,7 @@ endif()
 include(CheckFortranSourceCompiles)
 include(CheckFortranSourceRuns)
 check_fortran_source_compiles("error stop i; end" f18errorstop SRC_EXT f90)
+check_fortran_source_compiles("real, allocatable :: array(:, :, :, :, :, :, :, :, :, :); end" f03rank SRC_EXT f90)
 check_fortran_source_runs("use, intrinsic :: iso_fortran_env, only : real128; real(real128) :: x; x = x+1; end" f03real128)
 
 add_subdirectory(src)
@@ -1,5 +1,24 @@
 # Fortran Standard Library
 
+[![Actions Status](https://github.com/fortran-lang/stdlib/workflows/CI/badge.svg)](https://github.com/fortran-lang/stdlib/actions)
+[![Actions Status](https://github.com/fortran-lang/stdlib/workflows/CI_windows/badge.svg)](https://github.com/fortran-lang/stdlib/actions)
+
+
+## Goals and Motivation
+
+The Fortran Standard, as published by the ISO (https://wg5-fortran.org/), does
+not have a Standard Library. The goal of this project is to provide a community
+driven and agreed upon *de facto* "standard" library for Fortran, called a
+Fortran Standard Library (`stdlib`). We have a rigorous process how `stdlib` is
+developed as documented in our [Workflow](WORKFLOW.md). `stdlib` is both a
+specification and a reference implementation. We are cooperating with the
+Fortran Standards Committee (e.g., the effort
+[started](https://github.com/j3-fortran/fortran_proposals/issues/104) at the J3
+committee repository) and the plan is to continue working with the Committee in
+the future (such as in the step 5. in the [Workflow](WORKFLOW.md) document), so
+that if the Committee wants to standardize some feature already available in `stdlib`, it would
+base it on `stdlib`'s implementation.
+
 ## Scope
 
 The goal of the Fortran Standard Library is to achieve the following general scope:
@@ -16,19 +35,19 @@ The goal of the Fortran Standard Library is to achieve the following general sco
 
 ### Get the code
 
-```
+```sh
 git clone https://github.com/fortran-lang/stdlib
 cd stdlib
 ```
 
 ### Build with CMake
 
-```
-mkdir build
-cd build
-cmake ..
-make
-ctest
+```sh
+cmake -B build
+
+cmake --build build
+
+cmake --build build --target test
 ```
 
 ### Build with make
 
@@ -1,5 +1,57 @@
 # Workflow for the Fortran stdlib contributors
 
-This document will describe the workflow we'll follow when developing the Fortran stdlib.
-It's largely to be discussed and decided.
-For now, take a look at the [issues](https://github.com/fortran-lang/stdlib).
+This document describes our current workflow.
+
+We welcome everyone and anyone to participate and propose additions to stdlib.
+It is okay if you do not have experience for specification or implementation,
+but have an idea for stdlib. If the idea is popular among the community, more
+experienced contributors will help it through all 5 steps.
+
+
+1. **Idea**: You have an idea or a proposal. Open an
+   [issue](https://github.com/fortran-lang/stdlib/issues) to discuss it. This
+   is on the level of "is there interest in having image reader/writer
+   functions in stdlib?" The goal of this step is to find out if the community
+   is interested in having this functionality as part of stdlib.
+
+2. **API**: When there seems to be significant interest in the proposal (vast
+   majority of participants think it is a good idea), move on to discuss the
+   specific API. It's OK to propose the API off the bat if you already have an
+   idea for it. This step is exploratory and its goal is to find out what the
+   API should *look* and *feel* like.
+
+3. **Specification**: Discuss the API and iterate. When there is vast majority
+   approval for the API, move on to implement it and submit a PR. Small PRs are
+   always better than large.  It is OK to implement only a few functions of a
+   new module, and continue work on the others in a later PR. All new
+   functionality goes into an "experimental" namespace
+   (`stdlib_experimental_*.f90`).  As part of the PR, when submitting a new
+   public facing API, please provide the initial draft of the specification
+   document as well as the the initial reference implementation of this
+   specification.  The specification is a document that describes the API and
+   the functionality, so that anyone can use it to create an implementation
+   from scratch without looking at `stdlib`. The `stdlib` library then provides
+   the reference implementation.
+
+4. **Implementation** in experimental: When opening a PR, request reviews from
+   one or more people that are most relevant to it. These are likely to be
+   people involved in prior steps of the workflow. Other contributors (not
+   explicitly invited) are encouraged to provide reviews and suggestions as
+   well. Iterate until all (or most) participants are on the same page.
+   We can merge when there is vast majority approval of the PR.
+
+5. **Release**: Moving from experimental to release. The experimental
+   "namespace" contains new functionality together with its specification. In
+   order to move from experimental to release, the specification document must
+   be approved by the wide community and the standards committee (informally).
+   If that happens, it has now been blessed for broad use and we can move the
+   code into the main section of `stdlib`, and the particular specification
+   document becomes part of the Fortran Standard Library.
+
+
+Note: the general term "vast majority" above means at least 80%, but ultimately
+it is left to our best judgement to ensure that the community agrees that each
+PR and proposal was approved by "vast majority".
+
+You are welcome to propose changes to this workflow by opening an
+[issue](https://github.com/fortran-lang/stdlib/issues).
@@ -2,22 +2,40 @@ set(SRC
     stdlib_experimental_ascii.f90
     stdlib_experimental_io.f90
     stdlib_experimental_error.f90
-    stdlib_experimental_optval.f90
     stdlib_experimental_kinds.f90
+    stdlib_experimental_optval.f90
+    stdlib_experimental_system.F90
 )
 
 add_library(fortran_stdlib ${SRC})
 
+set(LIB_MOD_DIR ${CMAKE_CURRENT_BINARY_DIR}/mod_files/)
+set_target_properties(fortran_stdlib PROPERTIES
+    Fortran_MODULE_DIRECTORY ${LIB_MOD_DIR})
+target_include_directories(fortran_stdlib PUBLIC
+    $<BUILD_INTERFACE:${LIB_MOD_DIR}>
+    $<INSTALL_INTERFACE:include>
+)
+
 if(f18errorstop)
   target_sources(fortran_stdlib PRIVATE f18estop.f90)
 else()
   target_sources(fortran_stdlib PRIVATE f08estop.f90)
 endif()
 
+if(f03rank)
+  target_sources(fortran_stdlib PRIVATE f03_stdlib_experimental_stats.f90)
+  target_sources(fortran_stdlib PRIVATE f03_stdlib_experimental_stats_mean.f90)
+else()
+  target_sources(fortran_stdlib PRIVATE f90_stdlib_experimental_stats.f90)
+  target_sources(fortran_stdlib PRIVATE f90_stdlib_experimental_stats_mean.f90)
+endif()
+
 add_subdirectory(tests)
 
 install(TARGETS fortran_stdlib
         RUNTIME DESTINATION bin
         ARCHIVE DESTINATION lib
         LIBRARY DESTINATION lib
     )
+install(DIRECTORY ${LIB_MOD_DIR} DESTINATION include)
@@ -3,7 +3,9 @@ SRC = stdlib_experimental_ascii.f90 \
       stdlib_experimental_io.f90 \
       stdlib_experimental_optval.f90 \
       stdlib_experimental_kinds.f90 \
-      f18estop.f90
+      f18estop.f90 \
+      f03_stdlib_experimental_stats.f90 \
+      f03_stdlib_experimental_stats_mean.f90 
 
 LIB = libstdlib.a
 
@@ -34,3 +36,7 @@ stdlib_experimental_io.o: \
 	stdlib_experimental_optval.o \
 	stdlib_experimental_kinds.o
 stdlib_experimental_optval.o: stdlib_experimental_kinds.o
+f03_stdlib_experimental_stats_mean.o: \
+	stdlib_experimental_optval.o \
+	stdlib_experimental_kinds.o \
+	f03_stdlib_experimental_stats.o
@@ -218,9 +218,12 @@ subroutine qsavetxt(filename, d)
 ! call savetxt("log.txt", data)
 
 integer :: s, i
+character(len=14) :: format_string
+
+write(format_string, '(a1,i06,a7)') '(', size(d, 2), 'f40.34)'
 s = open(filename, "w")
 do i = 1, size(d, 1)
-    write(s, *) d(i, :)
+    write(s, format_string) d(i, :)
 end do
 close(s)
 end subroutine
@@ -332,15 +335,15 @@ integer function open(filename, mode, iostat) result(u)
 
 select case (mode_(3:3))
 case('t')
-    access_='sequential'
     form_='formatted'
-case('b', 's')
-    access_='stream'
+case('b')
     form_='unformatted'
 case default
     call error_stop("Unsupported mode: "//mode_(3:3))
 end select
 
+access_ = 'stream'
+
 if (present(iostat)) then
     open(newunit=u, file=filename, &
          action = action_, position = position_, status = status_, &
 
@@ -34,7 +34,7 @@ module stdlib_experimental_optval
 contains
 
 
-  pure function optval_sp(x, default) result(y)
+  pure elemental function optval_sp(x, default) result(y)
     real(sp), intent(in), optional :: x
     real(sp), intent(in) :: default
     real(sp) :: y
@@ -47,7 +47,7 @@ pure function optval_sp(x, default) result(y)
   end function optval_sp
 
 
-  pure function optval_dp(x, default) result(y)
+  pure elemental function optval_dp(x, default) result(y)
     real(dp), intent(in), optional :: x
     real(dp), intent(in) :: default
     real(dp) :: y
@@ -60,7 +60,7 @@ pure function optval_dp(x, default) result(y)
   end function optval_dp
 
 
-  pure function optval_qp(x, default) result(y)
+  pure elemental function optval_qp(x, default) result(y)
     real(qp), intent(in), optional :: x
     real(qp), intent(in) :: default
     real(qp) :: y
@@ -73,7 +73,7 @@ pure function optval_qp(x, default) result(y)
   end function optval_qp
 
 
-  pure function optval_int8(x, default) result(y)
+  pure elemental function optval_int8(x, default) result(y)
     integer(int8), intent(in), optional :: x
     integer(int8), intent(in) :: default
     integer(int8) :: y
@@ -86,7 +86,7 @@ pure function optval_int8(x, default) result(y)
   end function optval_int8
 
 
-  pure function optval_int16(x, default) result(y)
+  pure elemental function optval_int16(x, default) result(y)
     integer(int16), intent(in), optional :: x
     integer(int16), intent(in) :: default
     integer(int16) :: y
@@ -99,7 +99,7 @@ pure function optval_int16(x, default) result(y)
   end function optval_int16
 
 
-  pure function optval_int32(x, default) result(y)
+  pure elemental function optval_int32(x, default) result(y)
     integer(int32), intent(in), optional :: x
     integer(int32), intent(in) :: default
     integer(int32) :: y
@@ -112,7 +112,7 @@ pure function optval_int32(x, default) result(y)
   end function optval_int32
 
 
-  pure function optval_int64(x, default) result(y)
+  pure elemental function optval_int64(x, default) result(y)
     integer(int64), intent(in), optional :: x
     integer(int64), intent(in) :: default
     integer(int64) :: y
@@ -125,7 +125,7 @@ pure function optval_int64(x, default) result(y)
   end function optval_int64
 
 
-  pure function optval_logical(x, default) result(y)
+  pure elemental function optval_logical(x, default) result(y)
     logical, intent(in), optional :: x
     logical, intent(in) :: default
     logical :: y
 
@@ -0,0 +1,124 @@
+module stdlib_experimental_stats
+
+#:set VERSION90 = defined('VERSION90')
+#:set REALKINDS = ["sp", "dp", "qp"]
+#:set INTKINDS = ["int8", "int16", "int32", "int64"]
+#:set REALTYPES = ["real({})".format(k) for k in REALKINDS]
+#:set INTTYPES = ["integer({})".format(k) for k in INTKINDS]
+#:set iktr = list(zip(range(len(REALKINDS)), REALKINDS, REALTYPES))
+#:set ikti = list(zip(range(len(INTKINDS)), INTKINDS, INTTYPES))
+
+use stdlib_experimental_kinds, only: sp, dp, qp, &
+    int8, int16, int32, int64
+implicit none
+private
+! Public API
+public :: mean
+
+interface mean
+#:for i1, k1, t1 in iktr
+    module function mean_1_${k1}$_${k1}$(x) result(res)
+        ${t1}$, intent(in) :: x(:)
+        ${t1}$ :: res
+    end function mean_1_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_1_${k1}$_dp(x) result(res)
+        ${t1}$, intent(in) :: x(:)
+        real(dp) :: res
+    end function mean_1_${k1}$_dp
+#:endfor
+
+
+#:for i1, k1, t1 in iktr
+    module function mean_2_all_${k1}$_${k1}$(x) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        ${t1}$ :: res
+    end function mean_2_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_2_all_${k1}$_dp(x) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        real(dp) :: res
+    end function mean_2_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+    module function mean_2_${k1}$_${k1}$(x, dim) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        integer, intent(in) :: dim
+        ${t1}$ :: res(size(x)/size(x, dim))
+    end function mean_2_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+    module function mean_2_${k1}$_dp(x, dim) result(res)
+        ${t1}$, intent(in) :: x(:,:)
+        integer, intent(in) :: dim
+        real(dp) :: res(size(x)/size(x, dim))
+    end function mean_2_${k1}$_dp
+#:endfor
+
+
+#:def ranksuffix(rank)
+#{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
+#:enddef
+
+#:if VERSION90
+#:set ranks = range(3,8)
+#:else
+#:set ranks = range(3,16)
+#:endif
+
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    ${t1}$ :: res
+end function mean_${rank}$_all_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    real(dp) :: res
+end function mean_${rank}$_all_${k1}$_dp
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    ${t1}$ :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+end function mean_${rank}$_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    real(dp) :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+end function mean_${rank}$_${k1}$_dp
+#:endfor
+#:endfor
+
+end interface
+
+end module
@@ -0,0 +1,43 @@
+# Descriptive statistics
+
+## Implemented
+
+ * `mean`
+
+## MEAN - mean of array elements
+
+### Description:
+
+Returns the mean of all the elements of *array*, or of the elements of *array* along dimension *dim*.
+
+### Syntax:
+
+RESULT = mean(*array*)
+
+RESULT = mean(*array*, *dim*)
+
+### Arguments:
+
+*array*: Shall be an array of type INTEGER, or REAL.
+
+*dim* (optional): Shall be a scalar of type INTEGER with a value in the range from 1 to n, where n is the rank of *array*.
+
+### Return value:
+
+If *array* is of type REAL, the result is of the same type as array.
+If *array* is of type INTEGER, the result is of type as *double precision*.
+
+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.
+
+### Example:
+
+```fortran
+program test
+    use stdlib_experimental_stats, only: mean
+    implicit none
+    real :: x(1:6) = [ 1., 2., 3., 4., 5., 6. ]
+    print *, mean(x)                            !returns 21.
+    print *, mean( reshape(x, [ 2, 3 ] ))       !returns 21.
+    print *, mean( reshape(x, [ 2, 3 ] ), 1)    !returns [ 3., 7., 11. ]
+end program
+```
@@ -0,0 +1,293 @@
+submodule (stdlib_experimental_stats) stdlib_experimental_stats_mean
+
+#:set VERSION90 = defined('VERSION90')
+#:set REALKINDS = ["sp", "dp", "qp"]
+#:set INTKINDS = ["int8", "int16", "int32", "int64"]
+#:set REALTYPES = ["real({})".format(k) for k in REALKINDS]
+#:set INTTYPES = ["integer({})".format(k) for k in INTKINDS]
+#:set iktr = list(zip(range(len(REALKINDS)), REALKINDS, REALTYPES))
+#:set ikti = list(zip(range(len(INTKINDS)), INTKINDS, INTTYPES))
+
+use stdlib_experimental_error, only: error_stop
+implicit none
+
+contains
+
+#:for i1, k1, t1 in iktr
+module function mean_1_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:)
+    ${t1}$ :: res
+
+    integer :: i
+
+    res = 0.0_${k1}$
+    do i = 1, size(x)
+        res = res + x(i)
+    enddo
+    res = res / real(size(x), ${k1}$)
+
+end function mean_1_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_1_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:)
+    real(dp) :: res
+
+    integer :: i
+
+    res = 0.0_dp
+    do i = 1, size(x)
+        res = res + real(x(i), dp)
+    enddo
+    res = res / real(size(x), dp)
+
+end function mean_1_${k1}$_dp
+#:endfor
+
+
+#:for i1, k1, t1 in iktr
+module function mean_2_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    ${t1}$ :: res
+
+    integer :: i, i_
+
+    res = 0.0_${k1}$
+    do i_ = 1, size(x, 2)
+     do i = 1, size(x, 1)
+         res = res + x(i,i_)
+     enddo
+    enddo
+    res = res / real(size(x), ${k1}$)
+
+end function mean_2_all_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_2_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    real(dp) :: res
+
+    integer :: i, i_
+
+    res = 0.0_dp
+    do i_ = 1, size(x, 2)
+     do i = 1, size(x, 1)
+         res = res + real(x(i,i_), dp)
+     enddo
+    enddo
+    res = res / real(size(x), dp)
+
+end function mean_2_all_${k1}$_dp
+#:endfor
+
+#:for i1, k1, t1 in iktr
+module function mean_2_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    integer, intent(in) :: dim
+    ${t1}$ :: res(size(x)/size(x, dim))
+
+    integer :: i, i_
+
+    res = 0.0_${k1}$
+
+    select case(dim)
+     case(1)
+        do i_ = 1, size(x, 2)
+         do i = 1, size(x, 1)
+            res(i_) = res(i_) + x(i, i_)
+         end do
+        end do
+     case(2)
+        do i_ = 1, size(x, 2)
+         do i = 1, size(x, 1)
+            res(i) = res(i) + x(i, i_)
+         end do
+        end do
+     case default
+        call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+    res = res / real(size(x, dim), ${k1}$)
+
+end function mean_2_${k1}$_${k1}$
+#:endfor
+
+#:for i1, k1, t1 in ikti
+module function mean_2_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x(:,:)
+    integer, intent(in) :: dim
+    real(dp) :: res(size(x)/size(x, dim))
+
+    integer :: i, i_
+
+    res = 0.0_dp
+
+    select case(dim)
+     case(1)
+        do i_ = 1, size(x, 2)
+         do i = 1, size(x, 1)
+            res(i_) = res(i_) + real(x(i, i_), dp)
+         end do
+        end do
+     case(2)
+        do i_ = 1, size(x, 2)
+         do i = 1, size(x, 1)
+            res(i) = res(i) + real(x(i, i_), dp)
+         end do
+        end do
+     case default
+        call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+    res = res / real(size(x, dim), dp)
+
+end function mean_2_${k1}$_dp
+#:endfor
+
+
+#:def ranksuffix(rank)
+#{if rank > 0}#(${":" + ",:" * (rank - 1)}$)#{endif}#
+#:enddef
+
+#:def isuffix(rank)
+#:if rank > 0
+i#{for i in range(2,rank)}#${",i"+"_" * (i-1)}$#{endfor}#
+#:endif
+#:enddef
+
+#:def ressuffix(rank,dim)
+#:if rank > 0
+#{for i in range(1,dim-1)}#${"i" + "_" * (i-1)+","}$#{endfor}##{if dim-1 >0}#i${"_"*(dim-2) }$#{endif}##{if dim -1 >0 and dim <rank}#,#{endif}##{if dim <rank}#i${"_"*(dim) }$#{endif}##{for i in range(dim+1,rank)}#${",i"+"_" * (i)}$#{endfor}#
+#:endif
+#:enddef
+
+
+#:if VERSION90
+#:set ranks = range(3,8)
+#:else
+#:set ranks = range(3,16)
+#:endif
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_${k1}$(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    ${t1}$ :: res
+
+    integer :: ${isuffix(rank+1)}$
+
+    res = 0.0_${k1}$
+
+#:for fj in range(rank,0,-1)
+    ${" "* (rank - fj)}$do i${"_"* (fj-1)}$ = 1, size(x, ${fj}$)
+#:endfor
+    ${" "* (rank)}$res = res + x(${isuffix(rank+1)}$)
+#:for fj in range(rank,0,-1)
+    ${" "* (fj-1)}$end do
+#:endfor
+
+    res = res / real(size(x), ${k1}$)
+
+end function mean_${rank}$_all_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_all_${k1}$_dp(x) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    real(dp) :: res
+
+    integer :: ${isuffix(rank+1)}$
+
+    res = 0.0_dp
+
+#:for fj in range(rank,0,-1)
+    ${" "* (rank - fj)}$do i${"_"* (fj-1)}$ = 1, size(x, ${fj}$)
+#:endfor
+    ${" "* (rank)}$res = res + real(x(${isuffix(rank+1)}$), dp)
+#:for fj in range(rank,0,-1)
+    ${" "* (fj-1)}$end do
+#:endfor
+
+    res = res / real(size(x), dp)
+
+end function mean_${rank}$_all_${k1}$_dp
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in iktr
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_${k1}$(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    ${t1}$ :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+
+    integer :: ${isuffix(rank+1)}$
+
+    res = 0.0_${k1}$
+
+    select case(dim)
+#:for fi in range(1,rank+1)
+     case(${fi}$)
+#:for fj in range(rank,0,-1)
+       ${" "* (rank - fj)}$do i${"_"* (fj-1)}$ = 1, size(x, ${fj}$)
+#:endfor
+       ${" "* (rank)}$res(${ressuffix(rank,fi)}$) = res(${ressuffix(rank,fi)}$) + x(${isuffix(rank+1)}$)
+#:for fj in range(rank,0,-1)
+       ${" "* (fj-1)}$end do
+#:endfor
+#:endfor
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+    res = res / real(size(x, dim), ${k1}$)
+
+end function mean_${rank}$_${k1}$_${k1}$
+#:endfor
+#:endfor
+
+#:for i1, k1, t1 in ikti
+#:for rank in ranks
+module function mean_${rank}$_${k1}$_dp(x, dim) result(res)
+    ${t1}$, intent(in) :: x${ranksuffix(rank)}$
+    integer, intent(in) :: dim
+    real(dp) :: res( &
+#:for imerge in range(1,rank-1)
+                  merge(size(x,${imerge}$),size(x,${imerge + 1}$),mask = ${imerge}$ < dim ), &
+#:endfor
+                  merge(size(x,${rank-1}$),size(x,${rank}$),mask = ${rank-1}$ < dim ) )
+
+    integer :: ${isuffix(rank+1)}$
+
+    res = 0.0_dp
+
+    select case(dim)
+#:for fi in range(1,rank+1)
+     case(${fi}$)
+#:for fj in range(rank,0,-1)
+       ${" "* (rank - fj)}$do i${"_"* (fj-1)}$ = 1, size(x, ${fj}$)
+#:endfor
+       ${" "* (rank)}$res(${ressuffix(rank,fi)}$) = res(${ressuffix(rank,fi)}$) + real(x(${isuffix(rank+1)}$), dp)
+#:for fj in range(rank,0,-1)
+       ${" "* (fj-1)}$end do
+#:endfor
+#:endfor
+     case default
+       call error_stop("ERROR (mean): wrong dimension")
+    end select
+
+    res = res / real(size(x, dim), dp)
+
+end function mean_${rank}$_${k1}$_dp
+#:endfor
+#:endfor
+
+end submodule
@@ -0,0 +1,43 @@
+module stdlib_experimental_system
+use, intrinsic :: iso_c_binding, only : c_int, c_long
+implicit none
+private
+public :: sleep
+
+interface
+#ifdef _WIN32
+subroutine winsleep(dwMilliseconds) bind (C, name='Sleep')
+!! void Sleep(DWORD dwMilliseconds)
+!! https://docs.microsoft.com/en-us/windows/win32/api/synchapi/nf-synchapi-sleep
+import c_long
+integer(c_long), value, intent(in) :: dwMilliseconds
+end subroutine winsleep
+#else
+integer(c_int) function usleep(usec) bind (C)
+!! int usleep(useconds_t usec);
+!! https://linux.die.net/man/3/usleep
+import c_int
+integer(c_int), value, intent(in) :: usec
+end function usleep
+#endif
+end interface
+
+contains
+
+subroutine sleep(millisec)
+integer, intent(in) :: millisec
+integer(c_int) :: ierr
+
+#ifdef _WIN32
+!! PGI Windows, Ifort Windows, ....
+call winsleep(int(millisec, c_long))
+#else
+!! Linux, Unix, MacOS, MSYS2, ...
+ierr = usleep(int(millisec * 1000, c_int))
+if (ierr/=0) error stop 'problem with usleep() system call'
+#endif
+
+
+end subroutine sleep
+
+end module stdlib_experimental_system
@@ -9,6 +9,8 @@ endmacro(ADDTEST)
 add_subdirectory(ascii)
 add_subdirectory(io)
 add_subdirectory(optval)
+add_subdirectory(stats)
+add_subdirectory(system)
 
 ADDTEST(always_skip)
 set_tests_properties(always_skip PROPERTIES SKIP_RETURN_CODE 77)
 
@@ -4,13 +4,16 @@ all:
 	$(MAKE) -f Makefile.manual --directory=ascii
 	$(MAKE) -f Makefile.manual --directory=io
 	$(MAKE) -f Makefile.manual --directory=optval
+	$(MAKE) -f Makefile.manual --directory=stats
 
 test:
 	$(MAKE) -f Makefile.manual --directory=ascii test
 	$(MAKE) -f Makefile.manual --directory=io test
 	$(MAKE) -f Makefile.manual --directory=optval test
+	$(MAKE) -f Makefile.manual --directory=stats test
 
 clean:
 	$(MAKE) -f Makefile.manual --directory=ascii clean
 	$(MAKE) -f Makefile.manual --directory=io clean
 	$(MAKE) -f Makefile.manual --directory=optval clean
+	$(MAKE) -f Makefile.manual --directory=stats clean
@@ -20,16 +20,25 @@ program test_optval
 
   call test_optval_character
 
+
+  call test_optval_sp_arr
+  call test_optval_dp_arr
+  call test_optval_qp_arr
+
+  call test_optval_int8_arr
+  call test_optval_int16_arr
+  call test_optval_int32_arr
+  call test_optval_int64_arr
+
 contains
 
-  
   subroutine test_optval_sp
     print *, "test_optval_sp"
     call assert(foo_sp(1.0_sp) == 1.0_sp)
     call assert(foo_sp() == 2.0_sp)
   end subroutine test_optval_sp
 
-  
+
   function foo_sp(x) result(z)
     real(sp), intent(in), optional :: x
     real(sp) :: z
@@ -43,7 +52,7 @@ subroutine test_optval_dp
     call assert(foo_dp() == 2.0_dp)
   end subroutine test_optval_dp
 
-  
+
   function foo_dp(x) result(z)
     real(dp), intent(in), optional :: x
     real(dp) :: z
@@ -57,95 +66,207 @@ subroutine test_optval_qp
     call assert(foo_qp() == 2.0_qp)
   end subroutine test_optval_qp
 
-  
+
   function foo_qp(x) result(z)
     real(qp), intent(in), optional :: x
     real(qp) :: z
     z = optval(x, 2.0_qp)
   endfunction foo_qp
-  
-  
+
+
   subroutine test_optval_int8
     print *, "test_optval_int8"
     call assert(foo_int8(1_int8) == 1_int8)
     call assert(foo_int8() == 2_int8)
   end subroutine test_optval_int8
 
-  
+
   function foo_int8(x) result(z)
     integer(int8), intent(in), optional :: x
     integer(int8) :: z
     z = optval(x, 2_int8)
   endfunction foo_int8
-  
+
 
   subroutine test_optval_int16
     print *, "test_optval_int16"
     call assert(foo_int16(1_int16) == 1_int16)
     call assert(foo_int16() == 2_int16)
   end subroutine test_optval_int16
 
-  
+
   function foo_int16(x) result(z)
     integer(int16), intent(in), optional :: x
     integer(int16) :: z
     z = optval(x, 2_int16)
   endfunction foo_int16
 
-  
+
   subroutine test_optval_int32
     print *, "test_optval_int32"
     call assert(foo_int32(1_int32) == 1_int32)
     call assert(foo_int32() == 2_int32)
   end subroutine test_optval_int32
 
-  
+
   function foo_int32(x) result(z)
     integer(int32), intent(in), optional :: x
     integer(int32) :: z
     z = optval(x, 2_int32)
   endfunction foo_int32
 
-  
+
   subroutine test_optval_int64
     print *, "test_optval_int64"
     call assert(foo_int64(1_int64) == 1_int64)
     call assert(foo_int64() == 2_int64)
   end subroutine test_optval_int64
 
-  
+
   function foo_int64(x) result(z)
     integer(int64), intent(in), optional :: x
     integer(int64) :: z
     z = optval(x, 2_int64)
   endfunction foo_int64
-  
+
 
   subroutine test_optval_logical
     print *, "test_optval_logical"
     call assert(foo_logical(.true.))
     call assert(.not.foo_logical())
   end subroutine test_optval_logical
 
-  
+
   function foo_logical(x) result(z)
     logical, intent(in), optional :: x
     logical :: z
     z = optval(x, .false.)
   endfunction foo_logical
-  
+
 
   subroutine test_optval_character
     print *, "test_optval_character"
     call assert(foo_character("x") == "x")
     call assert(foo_character() == "y")
   end subroutine test_optval_character
 
-  
+
   function foo_character(x) result(z)
     character(len=*), intent(in), optional :: x
     character(len=:), allocatable :: z
     z = optval(x, "y")
   endfunction foo_character
-  
+
+
+  subroutine test_optval_sp_arr
+    print *, "test_optval_sp_arr"
+    call assert(all(foo_sp_arr([1.0_sp, -1.0_sp]) == [1.0_sp, -1.0_sp]))
+    call assert(all(foo_sp_arr() == [2.0_sp, -2.0_sp]))
+  end subroutine test_optval_sp_arr
+
+
+  function foo_sp_arr(x) result(z)
+    real(sp), dimension(2), intent(in), optional :: x
+    real(sp), dimension(2) :: z
+    z = optval(x, [2.0_sp, -2.0_sp])
+  end function foo_sp_arr
+
+
+  subroutine test_optval_dp_arr
+    print *, "test_optval_dp_arr"
+    call assert(all(foo_dp_arr([1.0_dp, -1.0_dp]) == [1.0_dp, -1.0_dp]))
+    call assert(all(foo_dp_arr() == [2.0_dp, -2.0_dp]))
+  end subroutine test_optval_dp_arr
+
+
+  function foo_dp_arr(x) result(z)
+    real(dp), dimension(2), intent(in), optional :: x
+    real(dp), dimension(2) :: z
+    z = optval(x, [2.0_dp, -2.0_dp])
+  end function foo_dp_arr
+
+
+  subroutine test_optval_qp_arr
+    print *, "test_optval_qp_arr"
+    call assert(all(foo_qp_arr([1.0_qp, -1.0_qp]) == [1.0_qp, -1.0_qp]))
+    call assert(all(foo_qp_arr() == [2.0_qp, -2.0_qp]))
+  end subroutine test_optval_qp_arr
+
+
+  function foo_qp_arr(x) result(z)
+    real(qp), dimension(2), intent(in), optional :: x
+    real(qp), dimension(2) :: z
+    z = optval(x, [2.0_qp, -2.0_qp])
+  end function foo_qp_arr
+
+
+  subroutine test_optval_int8_arr
+    print *, "test_optval_int8_arr"
+    call assert(all(foo_int8_arr([1_int8, -1_int8]) == [1_int8, -1_int8]))
+    call assert(all(foo_int8_arr() == [2_int8, -2_int8]))
+  end subroutine test_optval_int8_arr
+
+
+  function foo_int8_arr(x) result(z)
+    integer(int8), dimension(2), intent(in), optional :: x
+    integer(int8), dimension(2) :: z
+    z = optval(x, [2_int8, -2_int8])
+  end function foo_int8_arr
+
+
+  subroutine test_optval_int16_arr
+    print *, "test_optval_int16_arr"
+    call assert(all(foo_int16_arr([1_int16, -1_int16]) == [1_int16, -1_int16]))
+    call assert(all(foo_int16_arr() == [2_int16, -2_int16]))
+  end subroutine test_optval_int16_arr
+
+
+  function foo_int16_arr(x) result(z)
+    integer(int16), dimension(2), intent(in), optional :: x
+    integer(int16), dimension(2) :: z
+    z = optval(x, [2_int16, -2_int16])
+  end function foo_int16_arr
+
+
+  subroutine test_optval_int32_arr
+    print *, "test_optval_int32_arr"
+    call assert(all(foo_int32_arr([1_int32, -1_int32]) == [1_int32, -1_int32]))
+    call assert(all(foo_int32_arr() == [2_int32, -2_int32]))
+  end subroutine test_optval_int32_arr
+
+
+  function foo_int32_arr(x) result(z)
+    integer(int32), dimension(2), intent(in), optional :: x
+    integer(int32), dimension(2) :: z
+    z = optval(x, [2_int32, -2_int32])
+  end function foo_int32_arr
+
+
+  subroutine test_optval_int64_arr
+    print *, "test_optval_int64_arr"
+    call assert(all(foo_int64_arr([1_int64, -1_int64]) == [1_int64, -1_int64]))
+    call assert(all(foo_int64_arr() == [2_int64, -2_int64]))
+  end subroutine test_optval_int64_arr
+
+
+  function foo_int64_arr(x) result(z)
+    integer(int64), dimension(2), intent(in), optional :: x
+    integer(int64), dimension(2) :: z
+    z = optval(x, [2_int64, -2_int64])
+  end function foo_int64_arr
+
+
+  subroutine test_optval_logical_arr
+    print *, "test_optval_logical_arr"
+    call assert(all(foo_logical_arr()))
+    call assert(all(.not.foo_logical_arr()))
+  end subroutine test_optval_logical_arr
+
+
+  function foo_logical_arr(x) result(z)
+    logical, dimension(2), intent(in), optional :: x
+    logical, dimension(2) :: z
+    z = optval(x, [.false., .false.])
+  end function foo_logical_arr
+
 end program test_optval
@@ -0,0 +1,5 @@
+ADDTEST(mean)
+
+if(f03rank)
+ADDTEST(mean_f03)
+endif()
@@ -0,0 +1,3 @@
+PROGS_SRC = test_mean.f90
+
+include ../Makefile.manual.test.mk
@@ -0,0 +1,16 @@
+1.000000000000000021e-08 9.199998759392489944e+01
+1.024113254885563425e-08 9.199998731474968849e+01
+1.048233721895820948e-08 9.199998703587728244e+01
+1.072361403187881949e-08 9.199998675729767683e+01
+1.096496300919481796e-08 9.199998647900135040e+01
+1.120638417249036630e-08 9.199998620097916557e+01
+1.144787754335570897e-08 9.199998592322251056e+01
+1.168944314338753750e-08 9.199998564572304360e+01
+1.193108099418952317e-08 9.199998536847290609e+01
+1.217279111737088596e-08 9.199998509146449521e+01
+1.241457353454836993e-08 9.199998481469057765e+01
+1.265642826734443823e-08 9.199998453814424693e+01
+1.289835533738818635e-08 9.199998426181879552e+01
+1.314035476631514857e-08 9.199998398570787117e+01
+1.338242657576766519e-08 9.199998370980536322e+01
+1.362457078739434161e-08 9.199998343410533153e+01
@@ -0,0 +1,75 @@
+program test_mean
+use stdlib_experimental_error, only: assert
+use stdlib_experimental_kinds, only: sp, dp, int32, int64
+use stdlib_experimental_io, only: loadtxt
+use stdlib_experimental_stats, only: mean
+implicit none
+
+real(sp), allocatable :: s(:, :)
+real(dp), allocatable :: d(:, :)
+
+real(dp), allocatable :: d3(:, :, :)
+real(dp), allocatable :: d4(:, :, :, :)
+
+
+!sp
+call loadtxt("array3.dat", s)
+
+call assert( mean(s) - sum(s)/real(size(s), sp) == 0.0_sp)
+call assert( sum( abs( mean(s,1) - sum(s,1)/real(size(s,1), sp) )) == 0.0_sp)
+call assert( sum( abs( mean(s,2) - sum(s,2)/real(size(s,2), sp) )) == 0.0_sp)
+
+
+!dp
+call loadtxt("array3.dat", d)
+
+call assert( mean(d) - sum(d)/real(size(d), dp) == 0.0_dp)
+call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!int32
+call loadtxt("array3.dat", d)
+
+call assert( mean(int(d, int32)) - sum(real(int(d, int32),dp))/real(size(d), dp) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int32),1) - sum(real(int(d, int32),dp),1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int32),2) - sum(real(int(d, int32),dp),2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!int64
+call loadtxt("array3.dat", d)
+
+call assert( mean(int(d, int64)) - sum(real(int(d, int64),dp))/real(size(d), dp) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int64),1) - sum(real(int(d, int64),dp),1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum(abs( mean(int(d, int64),2) - sum(real(int(d, int64),dp),2)/real(size(d,2), dp) )) == 0.0_dp)
+
+
+!dp rank 3
+allocate(d3(size(d,1),size(d,2),3))
+d3(:,:,1)=d;
+d3(:,:,2)=d*1.5_dp;
+d3(:,:,3)=d*4._dp;
+
+call assert( mean(d3) - sum(d3)/real(size(d3), dp) == 0.0_dp)
+call assert( sum( abs( mean(d3,1) - sum(d3,1)/real(size(d3,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d3,2) - sum(d3,2)/real(size(d3,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d3,3) - sum(d3,3)/real(size(d3,3), dp) )) == 0.0_dp)
+
+
+!dp rank 4
+allocate(d4(size(d,1),size(d,2),3,9))
+d4 = 1.
+d4(:,:,1,1)=d;
+d4(:,:,2,1)=d*1.5_dp;
+d4(:,:,3,1)=d*4._dp;
+d4(:,:,3,9)=d*4._dp;
+
+call assert( mean(d4) - sum(d4)/real(size(d4), dp) == 0.0_dp)
+call assert( sum( abs( mean(d4,1) - sum(d4,1)/real(size(d4,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,2) - sum(d4,2)/real(size(d4,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,3) - sum(d4,3)/real(size(d4,3), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d4,4) - sum(d4,4)/real(size(d4,4), dp) )) == 0.0_dp)
+
+contains
+
+end program
@@ -0,0 +1,38 @@
+program test_mean
+use stdlib_experimental_error, only: assert
+use stdlib_experimental_kinds, only: sp, dp, int32, int64
+use stdlib_experimental_io, only: loadtxt
+use stdlib_experimental_stats, only: mean
+implicit none
+
+real(dp), allocatable :: d(:, :)
+real(dp), allocatable :: d8(:, :, :, :, :, :, :, :)
+
+
+!dp
+call loadtxt("array3.dat", d)
+
+call assert( mean(d) - sum(d)/real(size(d), dp) == 0.0_dp)
+call assert( sum( abs( mean(d,1) - sum(d,1)/real(size(d,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d,2) - sum(d,2)/real(size(d,2), dp) )) == 0.0_dp)
+
+!dp rank 8
+allocate(d8(size(d,1), size(d,2), 3, 4, 5, 6, 7, 8))
+d8(:, :, 1, 4, 5 ,6 ,7 ,8)=d;
+d8(:, :, 2, 4, 5 ,6 ,7 ,8)=d * 1.5_dp;
+d8(:, :, 3, 4, 5 ,6 ,7 ,8)=d * 4._dp;
+
+call assert( mean(d8) - sum(d8)/real(size(d8), dp) == 0.0_dp)
+
+call assert( sum( abs( mean(d8,1) - sum(d8,1)/real(size(d8,1), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,2) - sum(d8,2)/real(size(d8,2), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,3) - sum(d8,3)/real(size(d8,3), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,4) - sum(d8,4)/real(size(d8,4), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,5) - sum(d8,5)/real(size(d8,5), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,6) - sum(d8,6)/real(size(d8,6), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,7) - sum(d8,7)/real(size(d8,7), dp) )) == 0.0_dp)
+call assert( sum( abs( mean(d8,8) - sum(d8,8)/real(size(d8,8), dp) )) == 0.0_dp)
+
+contains
+
+end program
@@ -0,0 +1,5 @@
+add_executable(test_sleep test_sleep.f90)
+target_link_libraries(test_sleep fortran_stdlib)
+
+add_test(NAME Sleep COMMAND $<TARGET_FILE:test_sleep> 350)
+set_tests_properties(Sleep PROPERTIES TIMEOUT 1)
@@ -0,0 +1,33 @@
+program test_sleep
+use, intrinsic :: iso_fortran_env, only : int64, real64
+use stdlib_experimental_system, only : sleep
+
+implicit none
+
+integer :: ierr, millisec
+character(8) :: argv
+integer(int64) :: tic, toc, trate
+real(real64) :: t_ms
+
+call system_clock(count_rate=trate)
+
+millisec = 780
+call get_command_argument(1, argv, status=ierr)
+if (ierr==0) read(argv,*) millisec
+
+if (millisec<0) millisec=0
+
+call system_clock(count=tic)
+call sleep(millisec)
+call system_clock(count=toc)
+
+t_ms = (toc-tic) * 1000._real64 / trate
+
+if (millisec > 0) then
+  if (t_ms < 0.5 * millisec) error stop 'actual sleep time was too short'
+  if (t_ms > 2 * millisec) error stop 'actual sleep time was too long'
+endif
+
+print '(A,F8.3)', 'OK: test_sleep: slept for (ms): ',t_ms
+
+end program
Original file line number	Diff line number	Diff line change
`@@ -0,0 +1,3 @@`
	`1`	`+PROGS_SRC = test_mean.f90`
	`2`	`+`
	`3`	`+include ../Makefile.manual.test.mk`