diff --git a/doc/specs/stdlib_stats.md b/doc/specs/stdlib_stats.md
index b8a4f7f78..3f758a19d 100644
--- a/doc/specs/stdlib_stats.md
+++ b/doc/specs/stdlib_stats.md
@@ -204,9 +204,9 @@ and if `n` is an odd number, the median is:
median(array) = mean( array_sorted( floor( (n + 1) / 2.):floor( (n + 1) / 2.) + 1 ) )
```
-The current implementation is a quite naive implementation that relies on sorting
-the whole array, using the subroutine `[[stdlib_sorting(module):ord_sort(interface)]]`
-provided by the `[[stdlib_sorting(module)]]` module.
+The current implementation relies on a selection algorithm applied on a copy of
+the whole array, using the subroutine `[[stdlib_selection(module):select(interface)]]`
+provided by the `[[stdlib_selection(module)]]` module.
### Syntax
@@ -220,11 +220,11 @@ Generic subroutine
### Arguments
-`array`: Shall be an array of type `integer` or `real`.
+`array`: Shall be an array of type `integer` or `real`. It is an `intent(in)` argument.
-`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`. It is an `intent(in)` argument.
-`mask` (optional): Shall be of type `logical` and either 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`. It is an `intent(in)` argument.
### Return value
diff --git a/src/Makefile.manual b/src/Makefile.manual
index 2bfcb7091..3d02c4ec2 100644
--- a/src/Makefile.manual
+++ b/src/Makefile.manual
@@ -134,7 +134,7 @@ stdlib_stats_mean.o: \
stdlib_stats_median.o: \
stdlib_optval.o \
stdlib_kinds.o \
- stdlib_sorting.o \
+ stdlib_selection.o \
stdlib_stats.o
stdlib_stats_moment.o: \
stdlib_optval.o \
diff --git a/src/stdlib_stats_median.fypp b/src/stdlib_stats_median.fypp
index 5adb7b256..ef1e797ca 100644
--- a/src/stdlib_stats_median.fypp
+++ b/src/stdlib_stats_median.fypp
@@ -7,10 +7,7 @@ submodule (stdlib_stats) stdlib_stats_median
use, intrinsic:: ieee_arithmetic, only: ieee_value, ieee_quiet_nan, ieee_is_nan
use stdlib_error, only: error_stop
use stdlib_optval, only: optval
- ! Use "ord_sort" rather than "sort" because the former can be much faster for arrays
- ! that are already partly sorted. While it is slightly slower for random arrays,
- ! ord_sort seems a better overall choice.
- use stdlib_sorting, only: sort => ord_sort
+ use stdlib_selection, only: select
implicit none
contains
@@ -24,6 +21,7 @@ contains
real(${o1}$) :: res
integer(kind = int64) :: c, n
+ ${t1}$ :: val, val1
${t1}$, allocatable :: x_tmp(:)
if (.not.optval(mask, .true.) .or. size(x) == 0) then
@@ -43,16 +41,18 @@ contains
x_tmp = reshape(x, [n])
- call sort(x_tmp)
+ call select(x_tmp, c, val)
if (mod(n, 2_int64) == 0) then
+ val1 = minval(x_tmp(c+1:n)) !instead of call select(x_tmp, c+1, val1, left = c)
#:if t1[0] == 'r'
- res = sum(x_tmp(c:c+1)) / 2._${o1}$
+ res = (val + val1) / 2._${o1}$
#:else
- res = sum( real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+ res = (real(val, kind=${o1}$) + &
+ real(val1, kind=${o1}$)) / 2._${o1}$
#:endif
else
- res = x_tmp(c)
+ res = val
end if
end function ${name}$
@@ -74,6 +74,7 @@ contains
integer :: j${fj}$
#:endfor
#:endif
+ ${t1}$ :: val, val1
${t1}$, allocatable :: x_tmp(:)
if (.not.optval(mask, .true.) .or. size(x) == 0) then
@@ -107,17 +108,18 @@ contains
end if
#:endif
- call sort(x_tmp)
+ call select(x_tmp, c, val)
if (mod(n, 2) == 0) then
+ val1 = minval(x_tmp(c+1:n))
res${reduce_subvector('j', rank, fi)}$ = &
#:if t1[0] == 'r'
- sum(x_tmp(c:c+1)) / 2._${o1}$
+ (val + val1) / 2._${o1}$
#:else
- sum(real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
+ (real(val, kind=${o1}$) + real(val1, kind=${o1}$)) / 2._${o1}$
#:endif
else
- res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
+ res${reduce_subvector('j', rank, fi)}$ = val
end if
#:for fj in range(1, rank)
end do
@@ -141,6 +143,7 @@ contains
real(${o1}$) :: res
integer(kind = int64) :: c, n
+ ${t1}$ :: val, val1
${t1}$, allocatable :: x_tmp(:)
if (any(shape(x) .ne. shape(mask))) then
@@ -156,21 +159,26 @@ contains
x_tmp = pack(x, mask)
- call sort(x_tmp)
-
n = size(x_tmp, kind=int64)
- c = floor( (n + 1) / 2._${o1}$, kind=int64)
if (n == 0) then
res = ieee_value(1._${o1}$, ieee_quiet_nan)
- else if (mod(n, 2_int64) == 0) then
+ return
+ end if
+
+ c = floor( (n + 1) / 2._${o1}$, kind=int64)
+
+ call select(x_tmp, c, val)
+
+ if (mod(n, 2_int64) == 0) then
+ val1 = minval(x_tmp(c+1:n))
#:if t1[0] == 'r'
- res = sum(x_tmp(c:c+1)) / 2._${o1}$
+ res = (val + val1) / 2._${o1}$
#:else
- res = sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+ res = (real(val, kind=${o1}$) + real(val1, kind=${o1}$)) / 2._${o1}$
#:endif
else if (mod(n, 2_int64) == 1) then
- res = x_tmp(c)
+ res = val
end if
end function ${name}$
@@ -192,6 +200,7 @@ contains
integer :: j${fj}$
#:endfor
#:endif
+ ${t1}$ :: val, val1
${t1}$, allocatable :: x_tmp(:)
if (any(shape(x) .ne. shape(mask))) then
@@ -220,23 +229,28 @@ contains
end if
#:endif
- call sort(x_tmp)
-
n = size(x_tmp, kind=int64)
- c = floor( (n + 1) / 2._${o1}$, kind=int64 )
if (n == 0) then
res${reduce_subvector('j', rank, fi)}$ = &
ieee_value(1._${o1}$, ieee_quiet_nan)
- else if (mod(n, 2_int64) == 0) then
+ return
+ end if
+
+ c = floor( (n + 1) / 2._${o1}$, kind=int64 )
+
+ call select(x_tmp, c, val)
+
+ if (mod(n, 2_int64) == 0) then
+ val1 = minval(x_tmp(c+1:n))
res${reduce_subvector('j', rank, fi)}$ = &
#:if t1[0] == 'r'
- sum(x_tmp(c:c+1)) / 2._${o1}$
+ (val + val1) / 2._${o1}$
#:else
- sum(real(x_tmp(c:c+1), kind=${o1}$)) / 2._${o1}$
+ (real(val, kind=${o1}$) + real(val1, kind=${o1}$)) / 2._${o1}$
#:endif
else if (mod(n, 2_int64) == 1) then
- res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
+ res${reduce_subvector('j', rank, fi)}$ = val
end if
deallocate(x_tmp)
diff --git a/src/tests/stats/test_median.fypp b/src/tests/stats/test_median.fypp
index d257c7485..22d19785d 100644
--- a/src/tests/stats/test_median.fypp
+++ b/src/tests/stats/test_median.fypp
@@ -117,12 +117,6 @@ contains
call check(error, median(d2odd_${k1}$), 1._dp&
, 'median(d2odd_${k1}$): uncorrect answer'&
, thr = dptol)
- if (allocated(error)) return
-
- call check(error, median(d2odd_${k1}$), 1._dp&
- , 'median(d2odd_${k1}$): uncorrect answer'&
- , thr = dptol)
- if (allocated(error)) return
end subroutine