11#:include "common.fypp"
2+ #:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
23
34!! Licensing:
45!!
@@ -60,10 +61,10 @@ submodule(stdlib_sorting) stdlib_sorting_sort
6061contains
6162
6263
63- #:for k1 in INT_KINDS
64+ #:for k1, t1 in IR_KINDS_TYPES
6465
6566 pure module subroutine ${k1}$_sort( array )
66- ! `${k1}$_sort( array )` sorts the input `ARRAY` of type `INTEGER(${k1}$) `
67+ ! `${k1}$_sort( array )` sorts the input `ARRAY` of type `${t1}$ `
6768! using a hybrid sort based on the `introsort` of David Musser. As with
6869! `introsort`, `${k1}$_sort( array )` is an unstable hybrid comparison
6970! algorithm using `quicksort` for the main body of the sort tree,
@@ -73,11 +74,11 @@ contains
7374! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
7475! behavior is typically small compared to other sorting algorithms.
7576
76- integer(${k1}$) , intent(inout) :: array(0:)
77+ ${t1}$ , intent(inout) :: array(0:)
7778
7879 integer(int32) :: depth_limit
7980
80- depth_limit = 2 * int( floor( log( real( size( array, kind=int64 ), &
81+ depth_limit = 2 * int( floor( log( real( size( array, kind=int_size ), &
8182 kind=dp) ) / log(2.0_dp) ), &
8283 kind=int32 )
8384 call introsort(array, depth_limit)
@@ -89,7 +90,7 @@ contains
8990! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion
9091! of the `quicksort` is too slow as estimated from `DEPTH_LIMIT`,
9192! otherwise sorting is done by a `quicksort`.
92- integer(${k1}$) , intent(inout) :: array(0:)
93+ ${t1}$ , intent(inout) :: array(0:)
9394 integer(int32), intent(in) :: depth_limit
9495
9596 integer(int_size), parameter :: insert_size = 16_int_size
@@ -115,10 +116,10 @@ contains
115116
116117 pure subroutine partition( array, index )
117118! quicksort partition using median of three.
118- integer(${k1}$) , intent(inout) :: array(0:)
119+ ${t1}$ , intent(inout) :: array(0:)
119120 integer(int_size), intent(out) :: index
120121
121- integer(${k1}$) :: u, v, w, x, y
122+ ${t1}$ :: u, v, w, x, y
122123 integer(int_size) :: i, j
123124
124125! Determine median of three and exchange it with the end.
@@ -158,10 +159,10 @@ contains
158159
159160 pure subroutine insertion_sort( array )
160161! Bog standard insertion sort.
161- integer(${k1}$) , intent(inout) :: array(0:)
162+ ${t1}$ , intent(inout) :: array(0:)
162163
163164 integer(int_size) :: i, j
164- integer(${k1}$) :: key
165+ ${t1}$ :: key
165166
166167 do j=1_int_size, size(array, kind=int_size)-1
167168 key = array(j)
@@ -178,10 +179,10 @@ contains
178179
179180 pure subroutine heap_sort( array )
180181! A bog standard heap sort
181- integer(${k1}$) , intent(inout) :: array(0:)
182+ ${t1}$ , intent(inout) :: array(0:)
182183
183184 integer(int_size) :: i, heap_size
184- integer(${k1}$) :: y
185+ ${t1}$ :: y
185186
186187 heap_size = size( array, kind=int_size )
187188! Build the max heap
@@ -201,11 +202,11 @@ contains
201202
202203 pure recursive subroutine max_heapify( array, i, heap_size )
203204! Transform the array into a max heap
204- integer(${k1}$) , intent(inout) :: array(0:)
205+ ${t1}$ , intent(inout) :: array(0:)
205206 integer(int_size), intent(in) :: i, heap_size
206207
207208 integer(int_size) :: l, r, largest
208- integer(${k1}$) :: y
209+ ${t1}$ :: y
209210
210211 largest = i
211212 l = 2_int_size * i + 1_int_size
@@ -230,177 +231,6 @@ contains
230231#:endfor
231232
232233
233- #:for k1 in REAL_KINDS
234-
235- pure module subroutine ${k1}$_sort( array )
236-
237- ! `${k1}$_sort( array )` sorts the input `ARRAY` of type `REAL(${k1}$)`
238- ! using a hybrid sort based on the `introsort` of David Musser. As with
239- ! `introsort`, `${k1}$_sort( array )` is an unstable hybrid comparison
240- ! algorithm using `quicksort` for the main body of the sort tree,
241- ! supplemented by `insertion sort` for the outer brances, but if
242- ! `quicksort` is converging too slowly the algorithm resorts
243- ! to `heapsort`. The algorithm is of order O(N Ln(N)) for all inputs.
244- ! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
245- ! behavior is typically small compared to other sorting algorithms.
246-
247- real(${k1}$), intent(inout) :: array(0:)
248- integer(int32) :: depth_limit
249-
250- depth_limit = &
251- 2 * int( floor( log( real( size( array, kind=int_size ), &
252- kind=dp) ) / log(2.0_dp) ), &
253- kind=int32 )
254- call introsort(array, depth_limit)
255-
256- contains
257-
258- pure recursive subroutine introsort( array, depth_limit )
259- ! It devolves to `insertionsort` if the remaining number of elements
260- ! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion of
261- ! the quicksort is too slow as estimated from `DEPTH_LIMIT`, otherwise
262- ! sorting is done by a `quicksort`.
263- real(${k1}$), intent(inout) :: array(0:)
264- integer(int32), intent(in) :: depth_limit
265-
266- integer(int_size), parameter :: insert_size = 16_int_size
267- integer(int_size) :: index
268-
269- if ( size(array, kind=int_size) <= insert_size ) then
270- ! May be best at the end of SORT processing the whole array
271- ! See Musser, D.R., “Introspective Sorting and Selection
272- ! Algorithms,” Software—Practice and Experience, Vol. 27(8),
273- ! 983–993 (August 1997).
274-
275- call insertion_sort( array )
276-
277- else if ( depth_limit == 0 ) then
278- call heap_sort( array )
279-
280- else
281- call partition( array, index )
282- call introsort( array(0:index-1), depth_limit-1 )
283- call introsort( array(index+1:), depth_limit-1 )
284- end if
285-
286- end subroutine introsort
287-
288-
289- pure subroutine partition( array, index )
290- ! quicksort partition using median of three.
291- real(${k1}$), intent(inout) :: array(0:)
292- integer(int_size), intent(out) :: index
293-
294- real(${k1}$) :: u, v, w, x, y
295- integer(int_size) :: i, j
296-
297- ! Determine median of three and exchange it with the end.
298- u = array( 0 )
299- v = array( size(array, kind=int_size)/2-1 )
300- w = array( size(array, kind=int_size)-1 )
301- if ( (u > v) .neqv. (u > w) ) then
302- x = u
303- y = array(0)
304- array(0) = array( size( array, kind=int_size ) - 1 )
305- array( size( array, kind=int_size ) - 1 ) = y
306- else if ( (v < u) .neqv. (v < w) ) then
307- x = v
308- y = array(size( array, kind=int_size )/2-1)
309- array(size( array, kind=int_size )/2-1) = &
310- array(size( array, kind=int_size )-1)
311- array( size( array, kind=int_size )-1 ) = y
312- else
313- x = w
314- end if
315- ! Partition the array
316- i = -1_int_size
317- do j = 0_int_size, size(array, kind=int_size)-2
318- if ( array(j) <= x ) then
319- i = i + 1
320- y = array(i)
321- array(i) = array(j)
322- array(j) = y
323- end if
324- end do
325- y = array(i+1)
326- array(i+1) = array(size(array, kind=int_size)-1)
327- array(size(array, kind=int_size)-1) = y
328- index = i + 1
329-
330- end subroutine partition
331-
332- pure subroutine insertion_sort( array )
333- ! Bog standard insertion sort.
334- real(${k1}$), intent(inout) :: array(0:)
335-
336- integer(int_size) :: i, j
337- real(${k1}$) :: key
338-
339- do j=1_int_size, size(array, kind=int_size)-1
340- key = array(j)
341- i = j - 1
342- do while( i >= 0 )
343- if ( array(i) <= key ) exit
344- array(i+1) = array(i)
345- i = i - 1
346- end do
347- array(i+1) = key
348- end do
349-
350- end subroutine insertion_sort
351-
352- pure subroutine heap_sort( array )
353- ! A bog standard heap sort
354- real(${k1}$), intent(inout) :: array(0:)
355-
356- integer(int_size) :: i, heap_size
357- real(${k1}$) :: y
358-
359- heap_size = size( array, kind=int_size )
360- ! Build the max heap
361- do i = (heap_size-2)/2_int_size, 0_int_size, -1_int_size
362- call max_heapify( array, i, heap_size )
363- end do
364- do i = heap_size-1, 1_int_size, -1_int_size
365- ! Swap the first element with the current final element
366- y = array(0)
367- array(0) = array(i)
368- array(i) = y
369- ! Sift down using max_heapify
370- call max_heapify( array, 0_int_size, i )
371- end do
372-
373- end subroutine heap_sort
374-
375- pure recursive subroutine max_heapify( array, i, heap_size )
376- ! Transform the array into a max heap
377- real(${k1}$), intent(inout) :: array(0:)
378- integer(int_size), intent(in) :: i, heap_size
379-
380- integer(int_size) :: l, r, largest
381- real(${k1}$) :: y
382-
383- largest = i
384- l = 2_int_size * i + 1_int_size
385- r = l + 1_int_size
386- if ( l < heap_size ) then
387- if ( array(l) > array(largest) ) largest = l
388- end if
389- if ( r < heap_size ) then
390- if ( array(r) > array(largest) ) largest = r
391- end if
392- if ( largest /= i ) then
393- y = array(i)
394- array(i) = array(largest)
395- array(largest) = y
396- call max_heapify( array, largest, heap_size )
397- end if
398-
399- end subroutine max_heapify
400-
401- end subroutine ${k1}$_sort
402-
403- #:endfor
404234
405235
406236 pure module subroutine char_sort( array )
0 commit comments