1
1
#:include "common.fypp"
2
+ #:set IR_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES
2
3
3
4
!! Licensing:
4
5
!!
@@ -60,10 +61,10 @@ submodule(stdlib_sorting) stdlib_sorting_sort
60
61
contains
61
62
62
63
63
- #:for k1 in INT_KINDS
64
+ #:for k1, t1 in IR_KINDS_TYPES
64
65
65
66
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}$ `
67
68
! using a hybrid sort based on the `introsort` of David Musser. As with
68
69
! `introsort`, `${k1}$_sort( array )` is an unstable hybrid comparison
69
70
! algorithm using `quicksort` for the main body of the sort tree,
@@ -73,11 +74,11 @@ contains
73
74
! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
74
75
! behavior is typically small compared to other sorting algorithms.
75
76
76
- integer(${k1}$) , intent(inout) :: array(0:)
77
+ ${t1}$ , intent(inout) :: array(0:)
77
78
78
79
integer(int32) :: depth_limit
79
80
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 ), &
81
82
kind=dp) ) / log(2.0_dp) ), &
82
83
kind=int32 )
83
84
call introsort(array, depth_limit)
@@ -89,7 +90,7 @@ contains
89
90
! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion
90
91
! of the `quicksort` is too slow as estimated from `DEPTH_LIMIT`,
91
92
! otherwise sorting is done by a `quicksort`.
92
- integer(${k1}$) , intent(inout) :: array(0:)
93
+ ${t1}$ , intent(inout) :: array(0:)
93
94
integer(int32), intent(in) :: depth_limit
94
95
95
96
integer(int_size), parameter :: insert_size = 16_int_size
@@ -115,10 +116,10 @@ contains
115
116
116
117
pure subroutine partition( array, index )
117
118
! quicksort partition using median of three.
118
- integer(${k1}$) , intent(inout) :: array(0:)
119
+ ${t1}$ , intent(inout) :: array(0:)
119
120
integer(int_size), intent(out) :: index
120
121
121
- integer(${k1}$) :: u, v, w, x, y
122
+ ${t1}$ :: u, v, w, x, y
122
123
integer(int_size) :: i, j
123
124
124
125
! Determine median of three and exchange it with the end.
@@ -158,10 +159,10 @@ contains
158
159
159
160
pure subroutine insertion_sort( array )
160
161
! Bog standard insertion sort.
161
- integer(${k1}$) , intent(inout) :: array(0:)
162
+ ${t1}$ , intent(inout) :: array(0:)
162
163
163
164
integer(int_size) :: i, j
164
- integer(${k1}$) :: key
165
+ ${t1}$ :: key
165
166
166
167
do j=1_int_size, size(array, kind=int_size)-1
167
168
key = array(j)
@@ -178,10 +179,10 @@ contains
178
179
179
180
pure subroutine heap_sort( array )
180
181
! A bog standard heap sort
181
- integer(${k1}$) , intent(inout) :: array(0:)
182
+ ${t1}$ , intent(inout) :: array(0:)
182
183
183
184
integer(int_size) :: i, heap_size
184
- integer(${k1}$) :: y
185
+ ${t1}$ :: y
185
186
186
187
heap_size = size( array, kind=int_size )
187
188
! Build the max heap
@@ -201,11 +202,11 @@ contains
201
202
202
203
pure recursive subroutine max_heapify( array, i, heap_size )
203
204
! Transform the array into a max heap
204
- integer(${k1}$) , intent(inout) :: array(0:)
205
+ ${t1}$ , intent(inout) :: array(0:)
205
206
integer(int_size), intent(in) :: i, heap_size
206
207
207
208
integer(int_size) :: l, r, largest
208
- integer(${k1}$) :: y
209
+ ${t1}$ :: y
209
210
210
211
largest = i
211
212
l = 2_int_size * i + 1_int_size
@@ -230,177 +231,6 @@ contains
230
231
#:endfor
231
232
232
233
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
404
234
405
235
406
236
pure module subroutine char_sort( array )
0 commit comments