55
66module stdlib_selection
77!
8- ! This code was modified from a matlab implementation of "qselect " by Manolis
9- ! Lourakis , https://www.mathworks. com/matlabcentral/fileexchange/68947-qselect
10- ! Below is the licence of qselect
8+ ! This code was modified from the "Coretran" implementation "quickSelect " by
9+ ! Leon Folks , https://github. com/leonfoks/coretran/tree/master/src/sorting
10+ ! Below is the licence of "Coretran"
1111!
12- ! Copyright (c) 2018, Manolis Lourakis
12+ ! BSD 3-Clause License
13+ !
14+ ! Copyright (c) 2019, Leon Foks
1315! All rights reserved.
14- !
16+ !
1517! Redistribution and use in source and binary forms, with or without
1618! modification, are permitted provided that the following conditions are met:
17- !
18- ! * Redistributions of source code must retain the above copyright notice, this
19- ! list of conditions and the following disclaimer.
20- !
21- ! * Redistributions in binary form must reproduce the above copyright notice,
22- ! this list of conditions and the following disclaimer in the documentation
23- ! and/or other materials provided with the distribution
24- ! * Neither the name of Foundation for Research and Technology - Hellas nor the
25- ! names of its contributors may be used to endorse or promote products
26- ! derived from this software without specific prior written permission.
19+ !
20+ ! 1. Redistributions of source code must retain the above copyright notice, this
21+ ! list of conditions and the following disclaimer.
22+ !
23+ ! 2. Redistributions in binary form must reproduce the above copyright notice,
24+ ! this list of conditions and the following disclaimer in the documentation
25+ ! and/or other materials provided with the distribution.
26+ !
27+ ! 3. Neither the name of the copyright holder nor the names of its
28+ ! contributors may be used to endorse or promote products derived from
29+ ! this software without specific prior written permission.
30+ !
2731! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
2832! AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
2933! IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
30- ! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
34+ ! DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
3135! FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
3236! DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
3337! SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
3438! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
3539! OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
3640! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37- !
3841
3942use stdlib_kinds, only: int8, int16, int32, int64, sp, dp, qp
4043
@@ -70,13 +73,8 @@ contains
7073 subroutine ${name}$(a, k, kth_smallest, left, right)
7174 !! select - select the k-th smallest entry in a(:).
7275 !!
73- !! Implements Hoares Quickselect algorithm
74- !! (https://en.wikipedia.org/wiki/Quickselect)
75- !! with the median-of-3 pivot strategy.
76- !! Operates in-place, avoiding sorting and recursion.
77- !!
78- !! Modified from a translation of "qselect" by Manolis Lourakis
79- !! https://www.mathworks.com/matlabcentral/fileexchange/68947-qselect
76+ !! Partly derived from the "Coretran" implementation of
77+ !! quickSelect by Leon Folks, https://github.com/leonfoks/coretran
8078 !!
8179 ${arraytype}$, intent(inout) :: a(:)
8280 !! Array in which we seek the kth-smallest entry.
@@ -100,7 +98,7 @@ contains
10098 !! subroutine with different `k` (because of how `a(:)` becomes
10199 !! partially sorted, see documentation for `a(:)`).
102100
103- ${inttype}$ :: l, r, s, i, j, k_local
101+ ${inttype}$ :: l, r, mid, iPivot
104102 ${arraytype}$ :: pivot
105103 integer, parameter :: ip = ${intkind}$
106104
@@ -114,31 +112,20 @@ contains
114112 error stop "select must have 1 <= k <= size(a), and 1 <= left <= right <= size(a)";
115113 end if
116114
117- k_local = k - l + 1_ip
118-
119115 do while(.true.)
120- s = (l+r)/2_ip ! Deliberate integer division
121- if(a(s) < a(r)) call swap(a(s), a(r))
122- if(a(s) < a(l)) call swap(a(s), a(l))
123- if(a(r) < a(l)) call swap(a(r), a(l))
124- pivot = a(r) ! Median
125-
126- i = l
127- do j = l, r-1_ip
128- if(a(j) <= pivot) then
129- call swap(a(i), a(j))
130- i = i+1_ip
131- end if
132- end do
133- call swap(a(r), a(i))
134- s = i-l+1_ip
135- if(k_local < s) then
136- r = i-1_ip
137- else if(k_local > s) then
138- l=i+1_ip; k_local=k_local-s;
139- else
140- kth_smallest = a(i)
141- return
116+ mid = (l+r)/2_ip ! Deliberate integer division
117+
118+ call medianOf3(a, l, mid, r)
119+ call swap(a(l), a(mid))
120+ call partition(a, l, r, iPivot)
121+
122+ if (iPivot < k) then
123+ l = iPivot + 1
124+ elseif (iPivot > k) then
125+ r = iPivot - 1
126+ elseif (iPivot == k) then
127+ kth_smallest = a(k)
128+ return
142129 end if
143130 end do
144131
@@ -148,6 +135,41 @@ contains
148135 ${arraytype}$ :: tmp
149136 tmp = a; a = b; b = tmp
150137 end subroutine
138+
139+ subroutine medianOf3(a, left, mid, right)
140+ ${arraytype}$, intent(inout) :: a(:)
141+ ${inttype}$, intent(in) :: left, mid, right
142+ if(a(right) < a(left)) call swap(a(right), a(left))
143+ if(a(mid) < a(left)) call swap(a(mid) , a(left))
144+ if(a(right) < a(mid) ) call swap(a(mid) , a(right))
145+ end subroutine
146+
147+ subroutine partition(array,left,right,iPivot)
148+ ${arraytype}$, intent(inout) :: array(:)
149+ ${inttype}$, intent(in) :: left, right
150+ ${inttype}$, intent(out) :: iPivot
151+ ${inttype}$ :: lo,hi
152+ ${arraytype}$ :: pivot
153+
154+ pivot = array(left)
155+ lo = left
156+ hi=right
157+ do while (lo <= hi)
158+ do while (array(hi) > pivot)
159+ hi=hi-1
160+ end do
161+ do while (lo <= hi .and. array(lo) <= pivot)
162+ lo=lo+1
163+ end do
164+ if (lo <= hi) then
165+ call swap(array(lo),array(hi))
166+ lo=lo+1
167+ hi=hi-1
168+ end if
169+ end do
170+ call swap(array(left),array(hi))
171+ iPivot=hi
172+ end subroutine
151173 end subroutine
152174 #:endfor
153175 #:endfor
@@ -159,10 +181,8 @@ contains
159181 subroutine ${name}$(a, indx, k, kth_smallest, left, right)
160182 !! arg_select - find the index of the k-th smallest entry in `a(:)`
161183 !!
162- !! Implements Hoares Quickselect algorithm
163- !! https://en.wikipedia.org/wiki/Quickselect)
164- !! with the median-of-3 pivot strategy.
165- !! Operates in-place, avoiding sorting and recursion.
184+ !! Partly derived from the "Coretran" implementation of
185+ !! quickSelect by Leon Folks, https://github.com/leonfoks/coretran
166186 !!
167187 ${arraytype}$, intent(in) :: a(:)
168188 !! Array in which we seek the kth-smallest entry.
@@ -190,7 +210,7 @@ contains
190210 !! called the subroutine with a different `k` (due to the way that
191211 !! `indx(:)` becomes partially sorted, see documentation for `indx(:)`).
192212
193- ${inttype}$ :: l, r, s, i, j, k_local
213+ ${inttype}$ :: l, r, mid, iPivot
194214 ${arraytype}$ :: pivot
195215 integer, parameter :: ip = ${intkind}$
196216
@@ -208,31 +228,20 @@ contains
208228 error stop "arg_select must have 1 <= k <= size(a), and 1 <= left <= right <= size(a)";
209229 end if
210230
211- k_local = k - l + 1_ip
212-
213231 do while(.true.)
214- s = (l+r)/2_ip ! Deliberate integer division
215- if(a(indx(s)) < a(indx(r))) call swap(indx(s), indx(r))
216- if(a(indx(s)) < a(indx(l))) call swap(indx(s), indx(l))
217- if(a(indx(r)) < a(indx(l))) call swap(indx(r), indx(l))
218- pivot = a(indx(r)) ! Median
219-
220- i = l
221- do j = l, r-1_ip
222- if(a(indx(j)) <= pivot) then
223- call swap(indx(i), indx(j))
224- i = i+1_ip
225- end if
226- end do
227- call swap(indx(r), indx(i))
228- s = i-l+1_ip
229- if(k_local < s) then
230- r = i-1_ip
231- else if(k_local > s) then
232- l=i+1_ip; k_local=k_local-s;
233- else
234- kth_smallest = indx(i)
235- return
232+ mid = (l+r)/2_ip ! Deliberate integer division
233+
234+ call arg_medianOf3(a, indx, l, mid, r)
235+ call swap(indx(l), indx(mid))
236+ call arg_partition(a, indx, l, r, iPivot)
237+
238+ if (iPivot < k) then
239+ l = iPivot + 1
240+ elseif (iPivot > k) then
241+ r = iPivot - 1
242+ elseif (iPivot == k) then
243+ kth_smallest = indx(k)
244+ return
236245 end if
237246 end do
238247
@@ -242,6 +251,43 @@ contains
242251 ${inttype}$ :: tmp
243252 tmp = a; a = b; b = tmp
244253 end subroutine
254+
255+ subroutine arg_medianOf3(a, indx, left, mid, right)
256+ ${arraytype}$, intent(in) :: a(:)
257+ ${inttype}$, intent(inout) :: indx(:)
258+ ${inttype}$, intent(in) :: left, mid, right
259+ if(a(indx(right)) < a(indx(left))) call swap(indx(right), indx(left))
260+ if(a(indx(mid)) < a(indx(left))) call swap(indx(mid) , indx(left))
261+ if(a(indx(right)) < a(indx(mid)) ) call swap(indx(mid) , indx(right))
262+ end subroutine
263+
264+ subroutine arg_partition(array, indx, left,right,iPivot)
265+ ${arraytype}$, intent(in) :: array(:)
266+ ${inttype}$, intent(inout) :: indx(:)
267+ ${inttype}$, intent(in) :: left, right
268+ ${inttype}$, intent(out) :: iPivot
269+ ${inttype}$ :: lo,hi
270+ ${arraytype}$ :: pivot
271+
272+ pivot = array(indx(left))
273+ lo = left
274+ hi = right
275+ do while (lo <= hi)
276+ do while (array(indx(hi)) > pivot)
277+ hi=hi-1
278+ end do
279+ do while (lo <= hi .and. array(indx(lo)) <= pivot)
280+ lo=lo+1
281+ end do
282+ if (lo <= hi) then
283+ call swap(indx(lo),indx(hi))
284+ lo=lo+1
285+ hi=hi-1
286+ end if
287+ end do
288+ call swap(indx(left),indx(hi))
289+ iPivot=hi
290+ end subroutine
245291 end subroutine
246292 #:endfor
247293 #:endfor
0 commit comments