@@ -11,6 +11,7 @@ submodule (stdlib_stats) stdlib_stats_median
1111 ! that are already partly sorted. While it is slightly slower for random arrays,
1212 ! ord_sort seems a better overall choice.
1313 use stdlib_sorting, only: sort => ord_sort
14+ use stdlib_selection, only: select
1415 implicit none
1516
1617contains
@@ -24,6 +25,7 @@ contains
2425 real(${o1}$) :: res
2526
2627 integer(kind = int64) :: c, n
28+ ${t1}$ :: val, val1
2729 ${t1}$, allocatable :: x_tmp(:)
2830
2931 if (.not.optval(mask, .true.) .or. size(x) == 0) then
@@ -43,16 +45,18 @@ contains
4345
4446 x_tmp = reshape(x, [n])
4547
46- call sort (x_tmp)
48+ call select (x_tmp, c, val )
4749
4850 if (mod(n, 2_int64) == 0) then
51+ call select(x_tmp, c+1, val1)
4952 #:if t1[0] == 'r'
50- res = sum(x_tmp(c:c+1) ) / 2._${o1}$
53+ res = (val + val1 ) / 2._${o1}$
5154 #:else
52- res = sum( real(x_tmp(c:c+1), kind=${o1}$) ) / 2._${o1}$
55+ res = (real(val, kind=${o1}$) + &
56+ real(val1, kind=${o1}$)) / 2._${o1}$
5357 #:endif
5458 else
55- res = x_tmp(c)
59+ res = val
5660 end if
5761
5862 end function ${name}$
@@ -74,6 +78,7 @@ contains
7478 integer :: j${fj}$
7579 #:endfor
7680 #:endif
81+ ${t1}$ :: val, val1
7782 ${t1}$, allocatable :: x_tmp(:)
7883
7984 if (.not.optval(mask, .true.) .or. size(x) == 0) then
@@ -107,17 +112,18 @@ contains
107112 end if
108113 #:endif
109114
110- call sort (x_tmp)
115+ call select (x_tmp, c, val )
111116
112117 if (mod(n, 2) == 0) then
118+ call select(x_tmp, c+1, val1)
113119 res${reduce_subvector('j', rank, fi)}$ = &
114120 #:if t1[0] == 'r'
115- sum(x_tmp(c:c+1) ) / 2._${o1}$
121+ (val + val1 ) / 2._${o1}$
116122 #:else
117- sum (real(x_tmp(c:c+1) , kind=${o1}$) ) / 2._${o1}$
123+ (real(val , kind=${o1}$) + real(val1, kind=${o1}$) ) / 2._${o1}$
118124 #:endif
119125 else
120- res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
126+ res${reduce_subvector('j', rank, fi)}$ = val
121127 end if
122128 #:for fj in range(1, rank)
123129 end do
@@ -141,6 +147,7 @@ contains
141147 real(${o1}$) :: res
142148
143149 integer(kind = int64) :: c, n
150+ ${t1}$ :: val, val1
144151 ${t1}$, allocatable :: x_tmp(:)
145152
146153 if (any(shape(x) .ne. shape(mask))) then
@@ -156,21 +163,26 @@ contains
156163
157164 x_tmp = pack(x, mask)
158165
159- call sort(x_tmp)
160-
161166 n = size(x_tmp, kind=int64)
162- c = floor( (n + 1) / 2._${o1}$, kind=int64)
163167
164168 if (n == 0) then
165169 res = ieee_value(1._${o1}$, ieee_quiet_nan)
166- else if (mod(n, 2_int64) == 0) then
170+ return
171+ end if
172+
173+ c = floor( (n + 1) / 2._${o1}$, kind=int64)
174+
175+ call select(x_tmp, c, val)
176+
177+ if (mod(n, 2_int64) == 0) then
178+ call select(x_tmp, c+1, val1)
167179 #:if t1[0] == 'r'
168- res = sum(x_tmp(c:c+1) ) / 2._${o1}$
180+ res = (val + val1 ) / 2._${o1}$
169181 #:else
170- res = sum (real(x_tmp(c:c+1) , kind=${o1}$)) / 2._${o1}$
182+ res = (real(val, kind=${o1}$) + real(val1 , kind=${o1}$)) / 2._${o1}$
171183 #:endif
172184 else if (mod(n, 2_int64) == 1) then
173- res = x_tmp(c)
185+ res = val
174186 end if
175187
176188 end function ${name}$
@@ -192,6 +204,7 @@ contains
192204 integer :: j${fj}$
193205 #:endfor
194206 #:endif
207+ ${t1}$ :: val, val1
195208 ${t1}$, allocatable :: x_tmp(:)
196209
197210 if (any(shape(x) .ne. shape(mask))) then
@@ -220,23 +233,28 @@ contains
220233 end if
221234 #:endif
222235
223- call sort(x_tmp)
224-
225236 n = size(x_tmp, kind=int64)
226- c = floor( (n + 1) / 2._${o1}$, kind=int64 )
227237
228238 if (n == 0) then
229239 res${reduce_subvector('j', rank, fi)}$ = &
230240 ieee_value(1._${o1}$, ieee_quiet_nan)
231- else if (mod(n, 2_int64) == 0) then
241+ return
242+ end if
243+
244+ c = floor( (n + 1) / 2._${o1}$, kind=int64 )
245+
246+ call select(x_tmp, c, val)
247+
248+ if (mod(n, 2_int64) == 0) then
249+ call select(x_tmp, c+1, val1)
232250 res${reduce_subvector('j', rank, fi)}$ = &
233251 #:if t1[0] == 'r'
234- sum(x_tmp(c:c+1) ) / 2._${o1}$
252+ (val + val1 ) / 2._${o1}$
235253 #:else
236- sum (real(x_tmp(c:c+1) , kind=${o1}$)) / 2._${o1}$
254+ (real(val, kind=${o1}$) + real(val1 , kind=${o1}$)) / 2._${o1}$
237255 #:endif
238256 else if (mod(n, 2_int64) == 1) then
239- res${reduce_subvector('j', rank, fi)}$ = x_tmp(c)
257+ res${reduce_subvector('j', rank, fi)}$ = val
240258 end if
241259
242260 deallocate(x_tmp)
0 commit comments