Skip to content
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.

Commit 21e20f0

Browse files
jvdp1jvdp1
authoredSep 16, 2021
Generating sorting subroutines specific to character type with fypp (#475)
* sort_simp: remove subroutine specific to character, and are now generated by fypp
1 parent 608d5a2 commit 21e20f0

File tree

4 files changed

+104
-1052
lines changed

4 files changed

+104
-1052
lines changed
 

‎src/stdlib_sorting.fypp

Lines changed: 24 additions & 54 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,14 @@
11
#:include "common.fypp"
2-
#:set IRS_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + STRING_KINDS_TYPES
2+
3+
#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_KINDS))
4+
#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_KINDS))
5+
#:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_KINDS))
6+
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=len(array))"], ["char"]))
7+
8+
#! For better code reuse in fypp, make lists that contain the input types,
9+
#! with each having output types and a separate name prefix for subroutines
10+
#! This approach allows us to have the same code for all input types.
11+
#:set IRSC_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME
312

413
!! Licensing:
514
!!
@@ -353,30 +362,19 @@ module stdlib_sorting
353362
!! sorted data, having O(N) performance on uniformly non-increasing or
354363
!! non-decreasing data.
355364

356-
#:for k1, t1 in IRS_KINDS_TYPES
357-
module subroutine ${k1}$_ord_sort( array, work, reverse )
365+
#:for t1, t2, name1 in IRSC_TYPES_ALT_NAME
366+
module subroutine ${name1}$_ord_sort( array, work, reverse )
358367
!! Version: experimental
359368
!!
360-
!! `${k1}$_ord_sort( array )` sorts the input `ARRAY` of type `${t1}$`
369+
!! `${name1}$_ord_sort( array )` sorts the input `ARRAY` of type `${t1}$`
361370
!! using a hybrid sort based on the `'Rust" sort` algorithm found in `slice.rs`
362371
${t1}$, intent(inout) :: array(0:)
363-
${t1}$, intent(out), optional :: work(0:)
372+
${t2}$, intent(out), optional :: work(0:)
364373
logical, intent(in), optional :: reverse
365-
end subroutine ${k1}$_ord_sort
374+
end subroutine ${name1}$_ord_sort
366375

367376
#:endfor
368377

369-
module subroutine char_ord_sort( array, work, reverse )
370-
!! Version: experimental
371-
!!
372-
!! `char_ord_sort( array[, work, reverse] )` sorts the input `ARRAY` of type
373-
!! `CHARACTER(*)` using a hybrid sort based on the `'Rust" sort` algorithm
374-
!! found in `slice.rs`
375-
character(len=*), intent(inout) :: array(0:)
376-
character(len=len(array)), intent(out), optional :: work(0:)
377-
logical, intent(in), optional :: reverse
378-
end subroutine char_ord_sort
379-
380378
end interface ord_sort
381379

382380
interface sort
@@ -386,33 +384,21 @@ module stdlib_sorting
386384
!! on the `introsort` of David Musser.
387385
!! ([Specification](../page/specs/stdlib_sorting.html#sort-sorts-an-input-array))
388386

389-
#:for k1, t1 in IRS_KINDS_TYPES
390-
pure module subroutine ${k1}$_sort( array, reverse )
387+
#:for t1, t2, name1 in IRSC_TYPES_ALT_NAME
388+
pure module subroutine ${name1}$_sort( array, reverse )
391389
!! Version: experimental
392390
!!
393-
!! `${k1}$_sort( array[, reverse] )` sorts the input `ARRAY` of type `${t1}$`
391+
!! `${name1}$_sort( array[, reverse] )` sorts the input `ARRAY` of type `${t1}$`
394392
!! using a hybrid sort based on the `introsort` of David Musser.
395393
!! The algorithm is of order O(N Ln(N)) for all inputs.
396394
!! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
397395
!! behavior is small for random data compared to other sorting algorithms.
398396
${t1}$, intent(inout) :: array(0:)
399397
logical, intent(in), optional :: reverse
400-
end subroutine ${k1}$_sort
398+
end subroutine ${name1}$_sort
401399

402400
#:endfor
403401

404-
pure module subroutine char_sort( array, reverse )
405-
!! Version: experimental
406-
!!
407-
!! `char_sort( array[, reverse] )` sorts the input `ARRAY` of type
408-
!! `CHARACTER(*)` using a hybrid sort based on the `introsort` of David Musser.
409-
!! The algorithm is of order O(N Ln(N)) for all inputs.
410-
!! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
411-
!! behavior is small for random data compared to other sorting algorithms.
412-
character(len=*), intent(inout) :: array(0:)
413-
logical, intent(in), optional :: reverse
414-
end subroutine char_sort
415-
416402
end interface sort
417403

418404
interface sort_index
@@ -429,41 +415,25 @@ module stdlib_sorting
429415
!! non-decreasing sort, but if the optional argument `REVERSE` is present
430416
!! with a value of `.TRUE.` the indices correspond to a non-increasing sort.
431417

432-
#:for k1, t1 in IRS_KINDS_TYPES
433-
module subroutine ${k1}$_sort_index( array, index, work, iwork, &
418+
#:for t1, t2, name1 in IRSC_TYPES_ALT_NAME
419+
module subroutine ${name1}$_sort_index( array, index, work, iwork, &
434420
reverse )
435421
!! Version: experimental
436422
!!
437-
!! `${k1}$_sort_index( array, index[, work, iwork, reverse] )` sorts
423+
!! `${name1}$_sort_index( array, index[, work, iwork, reverse] )` sorts
438424
!! an input `ARRAY` of type `${t1}$`
439425
!! using a hybrid sort based on the `'Rust" sort` algorithm found in `slice.rs`
440426
!! and returns the sorted `ARRAY` and an array `INDEX of indices in the
441427
!! order that would sort the input `ARRAY` in the desired direction.
442428
${t1}$, intent(inout) :: array(0:)
443429
integer(int_size), intent(out) :: index(0:)
444-
${t1}$, intent(out), optional :: work(0:)
430+
${t2}$, intent(out), optional :: work(0:)
445431
integer(int_size), intent(out), optional :: iwork(0:)
446432
logical, intent(in), optional :: reverse
447-
end subroutine ${k1}$_sort_index
433+
end subroutine ${name1}$_sort_index
448434

449435
#:endfor
450436

451-
module subroutine char_sort_index( array, index, work, iwork, &
452-
reverse )
453-
!! Version: experimental
454-
!!
455-
!! `char_sort_index( array, index[, work, iwork, reverse] )` sorts
456-
!! an input `ARRAY` of type `CHARACTER(*)`
457-
!! using a hybrid sort based on the `'Rust" sort` algorithm found in `slice.rs`
458-
!! and returns the sorted `ARRAY` and an array `INDEX of indices in the
459-
!! order that would sort the input `ARRAY` in the desired direction.
460-
character(len=*), intent(inout) :: array(0:)
461-
integer(int_size), intent(out) :: index(0:)
462-
character(len=len(array)), intent(out), optional :: work(0:)
463-
integer(int_size), intent(out), optional :: iwork(0:)
464-
logical, intent(in), optional :: reverse
465-
end subroutine char_sort_index
466-
467437
end interface sort_index
468438

469439

‎src/stdlib_sorting_ord_sort.fypp

Lines changed: 32 additions & 378 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,13 @@
11
#:include "common.fypp"
2-
#:set IRS_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + STRING_KINDS_TYPES
2+
#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_TYPES, INT_KINDS))
3+
#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_TYPES, REAL_KINDS))
4+
#:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_TYPES, STRING_KINDS))
5+
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=:)"], ["character(len=len(array))"], ["char"]))
6+
7+
#! For better code reuse in fypp, make lists that contain the input types,
8+
#! with each having output types and a separate name prefix for subroutines
9+
#! This approach allows us to have the same code for all input types.
10+
#:set IRSC_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME
311

412
#:set SIGN_NAME = ["increase", "decrease"]
513
#:set SIGN_TYPE = [">", "<"]
@@ -61,10 +69,10 @@ submodule(stdlib_sorting) stdlib_sorting_ord_sort
6169
6270
contains
6371
64-
#:for k1, t1 in IRS_KINDS_TYPES
65-
module subroutine ${k1}$_ord_sort( array, work, reverse )
72+
#:for t1, t2, t3, name1 in IRSC_TYPES_ALT_NAME
73+
module subroutine ${name1}$_ord_sort( array, work, reverse )
6674
${t1}$, intent(inout) :: array(0:)
67-
${t1}$, intent(out), optional :: work(0:)
75+
${t3}$, intent(out), optional :: work(0:)
6876
logical, intent(in), optional :: reverse
6977
7078
logical :: reverse_
@@ -73,18 +81,18 @@ contains
7381
if(present(reverse)) reverse_ = reverse
7482
7583
if (reverse_) then
76-
call ${k1}$_decrease_ord_sort(array, work)
84+
call ${name1}$_decrease_ord_sort(array, work)
7785
else
78-
call ${k1}$_increase_ord_sort(array, work)
86+
call ${name1}$_increase_ord_sort(array, work)
7987
endif
8088
81-
end subroutine ${k1}$_ord_sort
89+
end subroutine ${name1}$_ord_sort
8290
#:endfor
8391
8492
#:for sname, signt, signoppt in SIGN_NAME_TYPE
85-
#:for k1, t1 in IRS_KINDS_TYPES
93+
#:for t1, t2, t3, name1 in IRSC_TYPES_ALT_NAME
8694
87-
subroutine ${k1}$_${sname}$_ord_sort( array, work )
95+
subroutine ${name1}$_${sname}$_ord_sort( array, work )
8896
! A translation to Fortran 2008, of the `"Rust" sort` algorithm found in
8997
! `slice.rs`
9098
! https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs#L2159
@@ -105,23 +113,28 @@ contains
105113
! original `listsort.txt`, and an optional `work` array to be used as
106114
! scratch memory.
107115
${t1}$, intent(inout) :: array(0:)
108-
${t1}$, intent(out), optional :: work(0:)
116+
${t3}$, intent(out), optional :: work(0:)
109117
110-
${t1}$, allocatable :: buf(:)
118+
${t2}$, allocatable :: buf(:)
111119
integer(int_size) :: array_size
112120
integer :: stat
113121
114122
array_size = size( array, kind=int_size )
115123
if ( present(work) ) then
116124
if ( size( work, kind=int_size) < array_size/2 ) then
117-
error stop "${k1}$_${sname}$_ord_sort: work array is too small."
125+
error stop "${name1}$_${sname}$_ord_sort: work array is too small."
118126
endif
119127
! Use the work array as scratch memory
120128
call merge_sort( array, work )
121129
else
122130
! Allocate a buffer to use as scratch memory.
131+
#:if t1[0:4] == "char"
132+
allocate( ${t3}$ :: buf(0:array_size/2-1), &
133+
stat=stat )
134+
#:else
123135
allocate( buf(0:array_size/2-1), stat=stat )
124-
if ( stat /= 0 ) error stop "${k1}$_${sname}$_ord_sort: Allocation of buffer failed."
136+
#:endif
137+
if ( stat /= 0 ) error stop "${name1}$_${sname}$_ord_sort: Allocation of buffer failed."
125138
call merge_sort( array, buf )
126139
end if
127140
@@ -153,7 +166,7 @@ contains
153166
${t1}$, intent(inout) :: array(0:)
154167
155168
integer(int_size) :: i, j
156-
${t1}$ :: key
169+
${t3}$ :: key
157170
158171
do j=1, size(array, kind=int_size)-1
159172
key = array(j)
@@ -229,7 +242,7 @@ contains
229242
230243
${t1}$, intent(inout) :: array(0:)
231244
232-
${t1}$ :: tmp
245+
${t3}$ :: tmp
233246
integer(int_size) :: i
234247
235248
tmp = array(0)
@@ -263,7 +276,7 @@ contains
263276
! worst-case.
264277
265278
${t1}$, intent(inout) :: array(0:)
266-
${t1}$, intent(inout) :: buf(0:)
279+
${t3}$, intent(inout) :: buf(0:)
267280
268281
integer(int_size) :: array_size, finish, min_run, r, r_count, &
269282
start
@@ -352,7 +365,7 @@ contains
352365
! must be long enough to hold the shorter of the two runs.
353366
${t1}$, intent(inout) :: array(0:)
354367
integer(int_size), intent(in) :: mid
355-
${t1}$, intent(inout) :: buf(0:)
368+
${t3}$, intent(inout) :: buf(0:)
356369

357370
integer(int_size) :: array_len, i, j, k
358371

@@ -408,7 +421,7 @@ contains
408421
${t1}$, intent(inout) :: array(0:)
409422

410423
integer(int_size) :: lo, hi
411-
${t1}$ :: temp
424+
${t3}$ :: temp
412425

413426
lo = 0
414427
hi = size( array, kind=int_size ) - 1
@@ -422,369 +435,10 @@ contains
422435

423436
end subroutine reverse_segment
424437

425-
end subroutine ${k1}$_${sname}$_ord_sort
438+
end subroutine ${name1}$_${sname}$_ord_sort
426439

427440
#:endfor
428441
#:endfor
429442

430-
module subroutine char_ord_sort( array, work, reverse )
431-
character(len=*), intent(inout) :: array(0:)
432-
character(len=len(array)), intent(out), optional :: work(0:)
433-
logical, intent(in), optional :: reverse
434-
435-
logical :: reverse_
436-
437-
reverse_ = .false.
438-
if(present(reverse)) reverse_ = reverse
439-
440-
if (reverse_) then
441-
call char_decrease_ord_sort(array, work)
442-
else
443-
call char_increase_ord_sort(array, work)
444-
endif
445-
446-
end subroutine char_ord_sort
447-
448-
449-
#:for sname, signt, signoppt in SIGN_NAME_TYPE
450-
subroutine char_${sname}$_ord_sort( array, work )
451-
! A translation to Fortran 2008, of the `"Rust" sort` algorithm found in
452-
! `slice.rs`
453-
! https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs#L2159
454-
! The Rust version in turn is a simplification of the Timsort algorithm
455-
! described in
456-
! https://svn.python.org/projects/python/trunk/Objects/listsort.txt, as
457-
! it drops both the use of 'galloping' to identify bounds of regions to be
458-
! sorted and the estimation of the optimal `run size`. However it remains
459-
! a hybrid sorting algorithm combining an iterative Merge sort controlled
460-
! by a stack of `RUNS` identified by regions of uniformly decreasing or
461-
! non-decreasing sequences that may be expanded to a minimum run size and
462-
! initially processed by an insertion sort.
463-
!
464-
! Note the Fortran implementation simplifies the logic as it only has to
465-
! deal with Fortran arrays of intrinsic types and not the full generality
466-
! of Rust's arrays and lists for arbitrary types. It also adds the
467-
! estimation of the optimal `run size` as suggested in Tim Peters'
468-
! original `listsort.txt`, and an optional `work` array to be used as
469-
! scratch memory.
470-
character(len=*), intent(inout) :: array(0:)
471-
character(len=len(array)), intent(out), optional :: work(0:)
472-
473-
character(len=:), allocatable :: buf(:)
474-
integer(int_size) :: array_size
475-
integer :: stat
476-
477-
if ( present(work) ) then
478-
! Use the work array as scratch memory
479-
call merge_sort( array, work )
480-
else
481-
! Allocate a buffer to use as scratch memory.
482-
array_size = size( array, kind=int_size )
483-
allocate( character(len=len(array)) :: buf(0:array_size/2-1), &
484-
stat=stat )
485-
if ( stat /= 0 ) error stop "${k1}$_${sname}$_ord_sort: Allocation of buffer failed."
486-
call merge_sort( array, buf )
487-
end if
488-
489-
contains
490-
491-
pure function calc_min_run( n ) result(min_run)
492-
!! Returns the minimum length of a run from 32-63 so that N/MIN_RUN is
493-
!! less than or equal to a power of two. See
494-
!! https://svn.python.org/projects/python/trunk/Objects/listsort.txt
495-
integer(int_size) :: min_run
496-
integer(int_size), intent(in) :: n
497-
498-
integer(int_size) :: num, r
499-
500-
num = n
501-
r = 0_int_size
502-
503-
do while( num >= 64 )
504-
r = ior( r, iand(num, 1_int_size) )
505-
num = ishft(num, -1_int_size)
506-
end do
507-
min_run = num + r
508-
509-
end function calc_min_run
510-
511-
512-
pure subroutine insertion_sort( array )
513-
! Sorts `ARRAY` using an insertion sort.
514-
character(len=*), intent(inout) :: array(0:)
515-
516-
integer(int_size) :: i, j
517-
character(len=len(array)) :: key
518-
519-
do j=1, size(array, kind=int_size)-1
520-
key = array(j)
521-
i = j - 1
522-
do while( i >= 0 )
523-
if ( array(i) ${signoppt}$= key ) exit
524-
array(i+1) = array(i)
525-
i = i - 1
526-
end do
527-
array(i+1) = key
528-
end do
529-
530-
end subroutine insertion_sort
531-
532-
533-
pure function collapse( runs ) result ( r )
534-
! Examine the stack of runs waiting to be merged, identifying adjacent runs
535-
! to be merged until the stack invariants are restablished:
536-
!
537-
! 1. len(-3) > len(-2) + len(-1)
538-
! 2. len(-2) > len(-1)
539-
integer(int_size) :: r
540-
type(run_type), intent(in), target :: runs(0:)
541-
542-
integer(int_size) :: n
543-
logical :: test
544-
545-
n = size(runs, kind=int_size)
546-
test = .false.
547-
if (n >= 2) then
548-
if ( runs( n-1 ) % base == 0 .or. &
549-
runs( n-2 ) % len <= runs(n-1) % len ) then
550-
test = .true.
551-
else if ( n >= 3 ) then ! X exists
552-
if ( runs(n-3) % len <= &
553-
runs(n-2) % len + runs(n-1) % len ) then
554-
test = .true.
555-
! |X| <= |Y| + |Z| => will need to merge due to rho1 or rho2
556-
else if( n >= 4 ) then
557-
if ( runs(n-4) % len <= &
558-
runs(n-3) % len + runs(n-2) % len ) then
559-
test = .true.
560-
! |W| <= |X| + |Y| => will need to merge due to rho1 or rho3
561-
end if
562-
end if
563-
end if
564-
end if
565-
if ( test ) then
566-
! By default merge Y & Z, rho2 or rho3
567-
if ( n >= 3 ) then
568-
if ( runs(n-3) % len < runs(n-1) % len ) then
569-
r = n - 3
570-
! |X| < |Z| => merge X & Y, rho1
571-
return
572-
end if
573-
end if
574-
r = n - 2
575-
! |Y| <= |Z| => merge Y & Z, rho4
576-
return
577-
else
578-
r = -1
579-
end if
580-
581-
end function collapse
582-
583-
584-
pure subroutine insert_head( array )
585-
! Inserts `array(0)` into the pre-sorted sequence `array(1:)` so that the
586-
! whole `array(0:)` becomes sorted, copying the first element into
587-
! a temporary variable, iterating until the right place for it is found.
588-
! copying every traversed element into the slot preceding it, and finally,
589-
! copying data from the temporary variable into the resulting hole.
590-
591-
character(len=*), intent(inout) :: array(0:)
592-
593-
character(len=len(array)) :: tmp
594-
integer(int_size) :: i
595-
596-
tmp = array(0)
597-
find_hole: do i=1, size(array, kind=int_size)-1
598-
if ( array(i) ${signt}$= tmp ) exit find_hole
599-
array(i-1) = array(i)
600-
end do find_hole
601-
array(i-1) = tmp
602-
603-
end subroutine insert_head
604-
605-
606-
subroutine merge_sort( array, buf )
607-
! The Rust merge sort borrows some (but not all) of the ideas from TimSort,
608-
! which is described in detail at
609-
! (http://svn.python.org/projects/python/trunk/Objects/listsort.txt).
610-
!
611-
! The algorithm identifies strictly descending and non-descending
612-
! subsequences, which are called natural runs. Where these runs are less
613-
! than a minimum run size they are padded by adding additional samples
614-
! using an insertion sort. The merge process is driven by a stack of
615-
! pending unmerged runs. Each newly found run is pushed onto the stack,
616-
! and then pairs of adjacentd runs are merged until these two invariants
617-
! are satisfied:
618-
!
619-
! 1. for every `i` in `1..size(runs)-1`: `runs(i - 1)%len > runs(i)%len`
620-
! 2. for every `i` in `2..size(runs)-1`: `runs(i - 2)%len >
621-
! runs(i - 1)%len + runs(i)%len`
622-
!
623-
! The invariants ensure that the total running time is `O(n log n)`
624-
! worst-case.
625-
626-
character(len=*), intent(inout) :: array(0:)
627-
character(len=len(array)), intent(inout) :: buf(0:)
628-
629-
integer(int_size) :: array_size, finish, min_run, r, r_count, &
630-
start
631-
type(run_type) :: runs(0:max_merge_stack-1), left, right
632-
633-
array_size = size(array, kind=int_size)
634-
635-
! Very short runs are extended using insertion sort to span at least
636-
! min_run elements. Slices of up to this length are sorted using insertion
637-
! sort.
638-
min_run = calc_min_run( array_size )
639-
640-
if ( array_size <= min_run ) then
641-
if ( array_size >= 2 ) call insertion_sort( array )
642-
return
643-
end if
644-
645-
! Following Rust sort, natural runs in `array` are identified by traversing
646-
! it backwards. By traversing it backward, merges more often go in the
647-
! opposite direction (forwards). According to developers of Rust sort,
648-
! merging forwards is slightly faster than merging backwards. Therefore
649-
! identifying runs by traversing backwards should improve performance.
650-
r_count = 0
651-
finish = array_size - 1
652-
do while ( finish >= 0 )
653-
! Find the next natural run, and reverse it if it's strictly descending.
654-
start = finish
655-
if ( start > 0 ) then
656-
start = start - 1
657-
if ( array(start+1) ${signoppt}$ array(start) ) then
658-
Descending: do while ( start > 0 )
659-
if ( array(start) ${signt}$= array(start-1) ) &
660-
exit Descending
661-
start = start - 1
662-
end do Descending
663-
call reverse_segment( array(start:finish) )
664-
else
665-
Ascending: do while( start > 0 )
666-
if ( array(start) ${signoppt}$ array(start-1) ) exit Ascending
667-
start = start - 1
668-
end do Ascending
669-
end if
670-
end if
671-
672-
! If the run is too short insert some more elements using an insertion sort.
673-
Insert: do while ( start > 0 )
674-
if ( finish - start >= min_run - 1 ) exit Insert
675-
start = start - 1
676-
call insert_head( array(start:finish) )
677-
end do Insert
678-
if ( start == 0 .and. finish == array_size - 1 ) return
679-
680-
runs(r_count) = run_type( base = start, &
681-
len = finish - start + 1 )
682-
finish = start-1
683-
r_count = r_count + 1
684-
685-
! Determine whether pairs of adjacent runs need to be merged to satisfy
686-
! the invariants, and, if so, merge them.
687-
Merge_loop: do
688-
r = collapse( runs(0:r_count - 1) )
689-
if ( r < 0 .or. r_count <= 1 ) exit Merge_loop
690-
left = runs( r + 1 )
691-
right = runs( r )
692-
call merge( array( left % base: &
693-
right % base + right % len - 1 ), &
694-
left % len, buf )
695-
696-
runs(r) = run_type( base = left % base, &
697-
len = left % len + right % len )
698-
if ( r == r_count - 3 ) runs(r+1) = runs(r+2)
699-
r_count = r_count - 1
700-
701-
end do Merge_loop
702-
end do
703-
if ( r_count /= 1 ) &
704-
error stop "MERGE_SORT completed without RUN COUNT == 1."
705-
706-
end subroutine merge_sort
707-
708-
709-
pure subroutine merge( array, mid, buf )
710-
! Merges the two non-decreasing runs `ARRAY(0:MID-1)` and `ARRAY(MID:)`
711-
! using `BUF` as temporary storage, and stores the merged runs into
712-
! `ARRAY(0:)`. `MID` must be > 0, and < `SIZE(ARRAY)-1`. Buffer `BUF`
713-
! must be long enough to hold the shorter of the two runs.
714-
character(len=*), intent(inout) :: array(0:)
715-
integer(int_size), intent(in) :: mid
716-
character(len=len(array)), intent(inout) :: buf(0:)
717-
718-
integer(int_size) :: array_len, i, j, k
719-
720-
array_len = size(array, kind=int_size)
721-
722-
! Merge first copies the shorter run into `buf`. Then, depending on which
723-
! run was shorter, it traces the copied run and the longer run forwards
724-
! (or backwards), comparing their next unprocessed elements and then
725-
! copying the lesser (or greater) one into `array`.
726-
727-
if ( mid <= array_len - mid ) then ! The left run is shorter.
728-
buf(0:mid-1) = array(0:mid-1)
729-
i = 0
730-
j = mid
731-
merge_lower: do k = 0, array_len-1
732-
if ( buf(i) ${signoppt}$= array(j) ) then
733-
array(k) = buf(i)
734-
i = i + 1
735-
if ( i >= mid ) exit merge_lower
736-
else
737-
array(k) = array(j)
738-
j = j + 1
739-
if ( j >= array_len ) then
740-
array(k+1:) = buf(i:mid-1)
741-
exit merge_lower
742-
end if
743-
end if
744-
end do merge_lower
745-
else ! The right run is shorter ! check that it is stable
746-
buf(0:array_len-mid-1) = array(mid:array_len-1)
747-
i = mid - 1
748-
j = array_len - mid -1
749-
merge_upper: do k = array_len-1, 0, -1
750-
if ( buf(j) ${signt}$= array(i) ) then
751-
array(k) = buf(j)
752-
j = j - 1
753-
if ( j < 0 ) exit merge_upper
754-
else
755-
array(k) = array(i)
756-
i = i - 1
757-
if ( i < 0 ) then
758-
array(0:k-1) = buf(0:j)
759-
exit merge_upper
760-
end if
761-
end if
762-
end do merge_upper
763-
end if
764-
end subroutine merge
765-
766-
767-
pure subroutine reverse_segment( array )
768-
! Reverse a segment of an array in place
769-
character(len=*), intent(inout) :: array(0:)
770-
771-
integer(int_size) :: lo, hi
772-
character(len=len(array)) :: temp
773-
774-
lo = 0
775-
hi = size( array, kind=int_size ) - 1
776-
do while( lo < hi )
777-
temp = array(lo)
778-
array(lo) = array(hi)
779-
array(hi) = temp
780-
lo = lo + 1
781-
hi = hi - 1
782-
end do
783-
784-
end subroutine reverse_segment
785-
786-
end subroutine char_${sname}$_ord_sort
787-
#:endfor
788-
789443
end submodule stdlib_sorting_ord_sort
790444

‎src/stdlib_sorting_sort.fypp

Lines changed: 23 additions & 201 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,13 @@
11
#:include "common.fypp"
2-
#:set IRS_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + STRING_KINDS_TYPES
2+
#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_KINDS))
3+
#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_KINDS))
4+
#:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_KINDS))
5+
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=len(array))"], ["char"]))
6+
7+
#! For better code reuse in fypp, make lists that contain the input types,
8+
#! with each having output types and a separate name prefix for subroutines
9+
#! This approach allows us to have the same code for all input types.
10+
#:set IRSC_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME
311

412
#:set SIGN_NAME = ["increase", "decrease"]
513
#:set SIGN_TYPE = [">", "<"]
@@ -65,8 +73,8 @@ submodule(stdlib_sorting) stdlib_sorting_sort
6573
6674
contains
6775
68-
#:for k1, t1 in IRS_KINDS_TYPES
69-
pure module subroutine ${k1}$_sort( array, reverse )
76+
#:for t1, t2, name1 in IRSC_TYPES_ALT_NAME
77+
pure module subroutine ${name1}$_sort( array, reverse )
7078
${t1}$, intent(inout) :: array(0:)
7179
logical, intent(in), optional :: reverse
7280
@@ -76,20 +84,20 @@ contains
7684
if(present(reverse)) reverse_ = reverse
7785
7886
if(reverse_)then
79-
call ${k1}$_decrease_sort(array)
87+
call ${name1}$_decrease_sort(array)
8088
else
81-
call ${k1}$_increase_sort(array)
89+
call ${name1}$_increase_sort(array)
8290
endif
83-
end subroutine ${k1}$_sort
91+
end subroutine ${name1}$_sort
8492
#:endfor
8593
8694
#:for sname, signt, signoppt in SIGN_NAME_TYPE
87-
#:for k1, t1 in IRS_KINDS_TYPES
95+
#:for t1, t2, name1 in IRSC_TYPES_ALT_NAME
8896
89-
pure subroutine ${k1}$_${sname}$_sort( array )
90-
! `${k1}$_${sname}$_sort( array )` sorts the input `ARRAY` of type `${t1}$`
97+
pure subroutine ${name1}$_${sname}$_sort( array )
98+
! `${name1}$_${sname}$_sort( array )` sorts the input `ARRAY` of type `${t1}$`
9199
! using a hybrid sort based on the `introsort` of David Musser. As with
92-
! `introsort`, `${k1}$_${sname}$_sort( array )` is an unstable hybrid comparison
100+
! `introsort`, `${name1}$_${sname}$_sort( array )` is an unstable hybrid comparison
93101
! algorithm using `quicksort` for the main body of the sort tree,
94102
! supplemented by `insertion sort` for the outer branches, but if
95103
! `quicksort` is converging too slowly the algorithm resorts
@@ -142,7 +150,7 @@ contains
142150
${t1}$, intent(inout) :: array(0:)
143151
integer(int_size), intent(out) :: index
144152
145-
${t1}$ :: u, v, w, x, y
153+
${t2}$ :: u, v, w, x, y
146154
integer(int_size) :: i, j
147155
148156
! Determine median of three and exchange it with the end.
@@ -185,7 +193,7 @@ contains
185193
${t1}$, intent(inout) :: array(0:)
186194
187195
integer(int_size) :: i, j
188-
${t1}$ :: key
196+
${t2}$ :: key
189197
190198
do j=1_int_size, size(array, kind=int_size)-1
191199
key = array(j)
@@ -205,7 +213,7 @@ contains
205213
${t1}$, intent(inout) :: array(0:)
206214
207215
integer(int_size) :: i, heap_size
208-
${t1}$ :: y
216+
${t2}$ :: y
209217
210218
heap_size = size( array, kind=int_size )
211219
! Build the max heap
@@ -229,7 +237,7 @@ contains
229237
integer(int_size), intent(in) :: i, heap_size
230238
231239
integer(int_size) :: l, r, largest
232-
${t1}$ :: y
240+
${t2}$ :: y
233241
234242
largest = i
235243
l = 2_int_size * i + 1_int_size
@@ -249,195 +257,9 @@ contains
249257
250258
end subroutine max_heapify
251259
252-
end subroutine ${k1}$_${sname}$_sort
260+
end subroutine ${name1}$_${sname}$_sort
253261
254262
#:endfor
255263
#:endfor
256264
257-
258-
pure module subroutine char_sort( array, reverse )
259-
character(len=*), intent(inout) :: array(0:)
260-
logical, intent(in), optional :: reverse
261-
262-
logical :: reverse_
263-
264-
reverse_ = .false.
265-
if(present(reverse)) reverse_ = reverse
266-
267-
if(reverse_)then
268-
call char_decrease_sort(array)
269-
else
270-
call char_increase_sort(array)
271-
endif
272-
end subroutine char_sort
273-
274-
275-
276-
#:for sname, signt, signoppt in SIGN_NAME_TYPE
277-
pure subroutine char_${sname}$_sort( array )
278-
! `char_${sname}$_sort( array )` sorts the input `ARRAY` of type `CHARACTER(*)`
279-
! using a hybrid sort based on the `introsort` of David Musser. As with
280-
! `introsort`, `char_${sname}$_sort( array )` is an unstable hybrid comparison
281-
! algorithm using `quicksort` for the main body of the sort tree,
282-
! supplemented by `insertion sort` for the outer branches, but if
283-
! `quicksort` is converging too slowly the algorithm resorts
284-
! to `heapsort`. The algorithm is of order O(N Ln(N)) for all inputs.
285-
! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
286-
! behavior is typically small compared to other sorting algorithms.
287-
288-
character(len=*), intent(inout) :: array(0:)
289-
290-
integer(int32) :: depth_limit
291-
292-
depth_limit = 2 * int( floor( log( real( size( array, kind=int_size ), &
293-
kind=dp) ) / log(2.0_dp) ), &
294-
kind=int32 )
295-
call introsort(array, depth_limit)
296-
297-
contains
298-
299-
pure recursive subroutine introsort( array, depth_limit )
300-
! It devolves to `insertionsort` if the remaining number of elements
301-
! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion
302-
! of the `quicksort` is too slow as estimated from `DEPTH_LIMIT`,
303-
! otherwise sorting is done by a `quicksort`.
304-
character(len=*), intent(inout) :: array(0:)
305-
integer(int32), intent(in) :: depth_limit
306-
307-
integer(int_size), parameter :: insert_size = 16_int_size
308-
integer(int_size) :: index
309-
310-
if ( size(array, kind=int_size) <= insert_size ) then
311-
! May be best at the end of SORT processing the whole array
312-
! See Musser, D.R., “Introspective Sorting and Selection
313-
! Algorithms,” Software—Practice and Experience, Vol. 27(8),
314-
! 983–993 (August 1997).
315-
316-
call insertion_sort( array )
317-
else if ( depth_limit == 0 ) then
318-
call heap_sort( array )
319-
else
320-
call partition( array, index )
321-
call introsort( array(0:index-1), depth_limit-1 )
322-
call introsort( array(index+1:), depth_limit-1 )
323-
end if
324-
325-
end subroutine introsort
326-
327-
328-
pure subroutine partition( array, index )
329-
! quicksort partition using median of three.
330-
character(len=*), intent(inout) :: array(0:)
331-
integer(int_size), intent(out) :: index
332-
333-
integer(int_size) :: i, j
334-
character(len=len(array)) :: u, v, w, x, y
335-
336-
! Determine median of three and exchange it with the end.
337-
u = array( 0 )
338-
v = array( size(array, kind=int_size)/2-1 )
339-
w = array( size(array, kind=int_size)-1 )
340-
if ( (u ${signt}$ v) .neqv. (u ${signt}$ w) ) then
341-
x = u
342-
y = array(0)
343-
array(0) = array( size( array, kind=int_size ) - 1 )
344-
array( size( array, kind=int_size ) - 1 ) = y
345-
else if ( (v ${signoppt}$ u) .neqv. (v ${signoppt}$ w) ) then
346-
x = v
347-
y = array(size( array, kind=int_size )/2-1)
348-
array( size( array, kind=int_size )/2-1 ) = &
349-
array( size( array, kind=int_size )-1 )
350-
array( size( array, kind=int_size )-1 ) = y
351-
else
352-
x = w
353-
end if
354-
! Partition the array.
355-
i = -1_int_size
356-
do j = 0_int_size, size(array, kind=int_size)-2
357-
if ( array(j) ${signoppt}$= x ) then
358-
i = i + 1
359-
y = array(i)
360-
array(i) = array(j)
361-
array(j) = y
362-
end if
363-
end do
364-
y = array(i+1)
365-
array(i+1) = array(size(array, kind=int_size)-1)
366-
array(size(array, kind=int_size)-1) = y
367-
index = i + 1
368-
369-
end subroutine partition
370-
371-
pure subroutine insertion_sort( array )
372-
! Bog standard insertion sort.
373-
character(len=*), intent(inout) :: array(0:)
374-
375-
integer(int_size) :: i, j
376-
character(len=len(array)) :: key
377-
378-
do j=1_int_size, size(array, kind=int_size)-1
379-
key = array(j)
380-
i = j - 1
381-
do while( i >= 0 )
382-
if ( array(i) ${signoppt}$= key ) exit
383-
array(i+1) = array(i)
384-
i = i - 1
385-
end do
386-
array(i+1) = key
387-
end do
388-
389-
end subroutine insertion_sort
390-
391-
pure subroutine heap_sort( array )
392-
! A bog standard heap sort
393-
character(len=*), intent(inout) :: array(0:)
394-
395-
integer(int_size) :: i, heap_size
396-
character(len=len(array)) :: y
397-
398-
heap_size = size( array, kind=int_size )
399-
! Build the max heap
400-
do i = (heap_size-2)/2_int_size, 0_int_size, -1_int_size
401-
call max_heapify( array, i, heap_size )
402-
end do
403-
do i = heap_size-1, 1_int_size, -1_int_size
404-
! Swap the first element with the current final element
405-
y = array(0)
406-
array(0) = array(i)
407-
array(i) = y
408-
! Sift down using max_heapify
409-
call max_heapify( array, 0_int_size, i )
410-
end do
411-
412-
end subroutine heap_sort
413-
414-
pure recursive subroutine max_heapify( array, i, heap_size )
415-
! Transform the array into a max heap
416-
character(len=*), intent(inout) :: array(0:)
417-
integer(int_size), intent(in) :: i, heap_size
418-
419-
integer(int_size) :: l, r, largest
420-
character(len=len(array)) :: y
421-
422-
largest = i
423-
l = 2_int_size * i + 1_int_size
424-
r = l + 1_int_size
425-
if ( l < heap_size ) then
426-
if ( array(l) ${signt}$ array(largest) ) largest = l
427-
end if
428-
if ( r < heap_size ) then
429-
if ( array(r) ${signt}$ array(largest) ) largest = r
430-
end if
431-
if ( largest /= i ) then
432-
y = array(i)
433-
array(i) = array(largest)
434-
array(largest) = y
435-
call max_heapify( array, largest, heap_size )
436-
end if
437-
438-
end subroutine max_heapify
439-
440-
end subroutine char_${sname}$_sort
441-
#:endfor
442-
443265
end submodule stdlib_sorting_sort

‎src/stdlib_sorting_sort_index.fypp

Lines changed: 25 additions & 419 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,13 @@
11
#:include "common.fypp"
2-
#:set IRS_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + STRING_KINDS_TYPES
2+
#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_TYPES, INT_KINDS))
3+
#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_TYPES, REAL_KINDS))
4+
#:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_TYPES, STRING_KINDS))
5+
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=:)"], ["character(len=len(array))"], ["char"]))
6+
7+
#! For better code reuse in fypp, make lists that contain the input types,
8+
#! with each having output types and a separate name prefix for subroutines
9+
#! This approach allows us to have the same code for all input types.
10+
#:set IRSC_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME
311

412
!! Licensing:
513
!!
@@ -56,10 +64,10 @@ submodule(stdlib_sorting) stdlib_sorting_sort_index
5664
5765
contains
5866
59-
#:for k1, t1 in IRS_KINDS_TYPES
67+
#:for t1, t2, t3, name1 in IRSC_TYPES_ALT_NAME
6068
61-
module subroutine ${k1}$_sort_index( array, index, work, iwork, reverse )
62-
! A modification of `${k1}$_ord_sort` to return an array of indices that
69+
module subroutine ${name1}$_sort_index( array, index, work, iwork, reverse )
70+
! A modification of `${name1}$_ord_sort` to return an array of indices that
6371
! would perform a stable sort of the `ARRAY` as input, and also sort `ARRAY`
6472
! as desired. The indices by default
6573
! correspond to a non-decreasing sort, but if the optional argument
@@ -85,12 +93,12 @@ contains
8593
8694
${t1}$, intent(inout) :: array(0:)
8795
integer(int_size), intent(out) :: index(0:)
88-
${t1}$, intent(out), optional :: work(0:)
96+
${t3}$, intent(out), optional :: work(0:)
8997
integer(int_size), intent(out), optional :: iwork(0:)
9098
logical, intent(in), optional :: reverse
9199
92100
integer(int_size) :: array_size, i, stat
93-
${t1}$, allocatable :: buf(:)
101+
${t2}$, allocatable :: buf(:)
94102
integer(int_size), allocatable :: ibuf(:)
95103
96104
array_size = size(array, kind=int_size)
@@ -121,7 +129,12 @@ contains
121129
call merge_sort( array, index, work, ibuf )
122130
end if
123131
else
132+
#:if t1[0:4] == "char"
133+
allocate( ${t3}$ :: buf(0:array_size/2-1), &
134+
stat=stat )
135+
#:else
124136
allocate( buf(0:array_size/2-1), stat=stat )
137+
#:endif
125138
if ( stat /= 0 ) error stop "Allocation of array buffer failed."
126139
if ( present(iwork) ) then
127140
if ( size(iwork, kind=int_size) < array_size/2 ) then
@@ -171,7 +184,7 @@ contains
171184
integer(int_size), intent(inout) :: index(0:)
172185
173186
integer(int_size) :: i, j, key_index
174-
${t1}$ :: key
187+
${t3}$ :: key
175188
176189
do j=1, size(array, kind=int_size)-1
177190
key = array(j)
@@ -254,7 +267,7 @@ contains
254267
${t1}$, intent(inout) :: array(0:)
255268
integer(int_size), intent(inout) :: index(0:)
256269
257-
${t1}$ :: tmp
270+
${t3}$ :: tmp
258271
integer(int_size) :: i, tmp_index
259272
260273
tmp = array(0)
@@ -293,7 +306,7 @@ contains
293306
294307
${t1}$, intent(inout) :: array(0:)
295308
integer(int_size), intent(inout) :: index(0:)
296-
${t1}$, intent(inout) :: buf(0:)
309+
${t3}$, intent(inout) :: buf(0:)
297310
integer(int_size), intent(inout) :: ibuf(0:)
298311
299312
integer(int_size) :: array_size, finish, min_run, r, r_count, &
@@ -386,7 +399,7 @@ contains
386399
! must be long enough to hold the shorter of the two runs.
387400
${t1}$, intent(inout) :: array(0:)
388401
integer(int_size), intent(in) :: mid
389-
${t1}$, intent(inout) :: buf(0:)
402+
${t3}$, intent(inout) :: buf(0:)
390403
integer(int_size), intent(inout) :: index(0:)
391404
integer(int_size), intent(inout) :: ibuf(0:)
392405

@@ -453,7 +466,7 @@ contains
453466
integer(int_size), intent(inout) :: index(0:)
454467

455468
integer(int_size) :: itemp, lo, hi
456-
${t1}$ :: temp
469+
${t3}$ :: temp
457470

458471
lo = 0
459472
hi = size( array, kind=int_size ) - 1
@@ -470,415 +483,8 @@ contains
470483

471484
end subroutine reverse_segment
472485

473-
end subroutine ${k1}$_sort_index
486+
end subroutine ${name1}$_sort_index
474487

475488
#:endfor
476489

477-
478-
module subroutine char_sort_index( array, index, work, iwork, reverse )
479-
! A modification of `char_ord_sort` to return an array of indices that
480-
! would perform a stable sort of the `ARRAY` as input, and also sort `ARRAY`
481-
! as desired. The indices by default
482-
! correspond to a non-decreasing sort, but if the optional argument
483-
! `REVERSE` is present with a value of `.TRUE.` the indices correspond to
484-
! a non-increasing sort. The logic of the determination of indexing largely
485-
! follows the `"Rust" sort` found in `slice.rs`:
486-
! https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs#L2159
487-
! The Rust version is a simplification of the Timsort algorithm described
488-
! in https://svn.python.org/projects/python/trunk/Objects/listsort.txt, as
489-
! it drops both the use of 'galloping' to identify bounds of regions to be
490-
! sorted and the estimation of the optimal `run size`. However it remains
491-
! a hybrid sorting algorithm combining an iterative Merge sort controlled
492-
! by a stack of `RUNS` identified by regions of uniformly decreasing or
493-
! non-decreasing sequences that may be expanded to a minimum run size, with
494-
! an insertion sort.
495-
!
496-
! Note the Fortran implementation simplifies the logic as it only has to
497-
! deal with Fortran arrays of intrinsic types and not the full generality
498-
! of Rust's arrays and lists for arbitrary types. It also adds the
499-
! estimation of the optimal `run size` as suggested in Tim Peters'
500-
! original listsort.txt, and the optional `work` and `iwork` arraya to be
501-
! used as scratch memory.
502-
503-
character(len=*), intent(inout) :: array(0:)
504-
integer(int_size), intent(out) :: index(0:)
505-
character(len=len(array)), intent(out), optional :: work(0:)
506-
integer(int_size), intent(out), optional :: iwork(0:)
507-
logical, intent(in), optional :: reverse
508-
509-
integer(int_size) :: array_size, i, stat
510-
character(len=:), allocatable :: buf(:)
511-
integer(int_size), allocatable :: ibuf(:)
512-
513-
array_size = size(array, kind=int_size)
514-
515-
do i = 0, array_size-1
516-
index(i) = i+1
517-
end do
518-
519-
if ( present(reverse) ) then
520-
if ( reverse ) then
521-
call reverse_segment( array, index )
522-
end if
523-
end if
524-
525-
! If necessary allocate buffers to serve as scratch memory.
526-
if ( present(work) ) then
527-
if ( present(iwork) ) then
528-
call merge_sort( array, index, work, iwork )
529-
else
530-
allocate( ibuf(0:array_size/2-1), stat=stat )
531-
if ( stat /= 0 ) error stop "Allocation of index buffer failed."
532-
call merge_sort( array, index, work, ibuf )
533-
end if
534-
else
535-
allocate( character(len=len(array)) :: buf(0:array_size/2-1), &
536-
stat=stat )
537-
if ( stat /= 0 ) error stop "Allocation of array buffer failed."
538-
if ( present(iwork) ) then
539-
call merge_sort( array, index, buf, iwork )
540-
else
541-
allocate( ibuf(0:array_size/2-1), stat=stat )
542-
if ( stat /= 0 ) error stop "Allocation of index buffer failed."
543-
call merge_sort( array, index, buf, ibuf )
544-
end if
545-
end if
546-
547-
if ( present(reverse) ) then
548-
if ( reverse ) then
549-
call reverse_segment( array, index )
550-
end if
551-
end if
552-
553-
contains
554-
555-
pure function calc_min_run( n ) result(min_run)
556-
!! Returns the minimum length of a run from 32-63 so that N/MIN_RUN is
557-
!! less than or equal to a power of two. See
558-
!! https://svn.python.org/projects/python/trunk/Objects/listsort.txt
559-
integer(int_size) :: min_run
560-
integer(int_size), intent(in) :: n
561-
562-
integer(int_size) :: num, r
563-
564-
num = n
565-
r = 0_int_size
566-
567-
do while( num >= 64 )
568-
r = ior( r, iand(num, 1_int_size) )
569-
num = ishft(num, -1_int_size)
570-
end do
571-
min_run = num + r
572-
573-
end function calc_min_run
574-
575-
576-
pure subroutine insertion_sort( array, index )
577-
! Sorts `ARRAY` using an insertion sort, while maintaining consistency in
578-
! location of the indices in `INDEX` to the elements of `ARRAY`.
579-
character(len=*), intent(inout) :: array(0:)
580-
integer(int_size), intent(inout) :: index(0:)
581-
582-
integer(int_size) :: i, j, key_index
583-
character(len=len(array)) :: key
584-
585-
do j=1, size(array, kind=int_size)-1
586-
key = array(j)
587-
key_index = index(j)
588-
i = j - 1
589-
do while( i >= 0 )
590-
if ( array(i) <= key ) exit
591-
array(i+1) = array(i)
592-
index(i+1) = index(i)
593-
i = i - 1
594-
end do
595-
array(i+1) = key
596-
index(i+1) = key_index
597-
end do
598-
599-
end subroutine insertion_sort
600-
601-
602-
pure function collapse( runs ) result ( r )
603-
! Examine the stack of runs waiting to be merged, identifying adjacent runs
604-
! to be merged until the stack invariants are restablished:
605-
!
606-
! 1. len(-3) > len(-2) + len(-1)
607-
! 2. len(-2) > len(-1)
608-
609-
integer(int_size) :: r
610-
type(run_type), intent(in), target :: runs(0:)
611-
612-
integer(int_size) :: n
613-
logical :: test
614-
615-
n = size(runs, kind=int_size)
616-
test = .false.
617-
if (n >= 2) then
618-
if ( runs( n-1 ) % base == 0 .or. &
619-
runs( n-2 ) % len <= runs(n-1) % len ) then
620-
test = .true.
621-
else if ( n >= 3 ) then ! X exists
622-
if ( runs(n-3) % len <= &
623-
runs(n-2) % len + runs(n-1) % len ) then
624-
test = .true.
625-
! |X| <= |Y| + |Z| => will need to merge due to rho1 or rho2
626-
else if( n >= 4 ) then
627-
if ( runs(n-4) % len <= &
628-
runs(n-3) % len + runs(n-2) % len ) then
629-
test = .true.
630-
! |W| <= |X| + |Y| => will need to merge due to rho1 or rho3
631-
end if
632-
end if
633-
end if
634-
end if
635-
if ( test ) then
636-
! By default merge Y & Z, rho2 or rho3
637-
if ( n >= 3 ) then
638-
if ( runs(n-3) % len < runs(n-1) % len ) then
639-
r = n - 3
640-
! |X| < |Z| => merge X & Y, rho1
641-
return
642-
end if
643-
end if
644-
r = n - 2
645-
! |Y| <= |Z| => merge Y & Z, rho4
646-
return
647-
else
648-
r = -1
649-
end if
650-
651-
end function collapse
652-
653-
654-
pure subroutine insert_head( array, index )
655-
! Inserts `array(0)` into the pre-sorted sequence `array(1:)` so that the
656-
! whole `array(0:)` becomes sorted, copying the first element into
657-
! a temporary variable, iterating until the right place for it is found.
658-
! copying every traversed element into the slot preceding it, and finally,
659-
! copying data from the temporary variable into the resulting hole.
660-
! Consistency of the indices in `index` with the elements of `array`
661-
! are maintained.
662-
663-
character(len=*), intent(inout) :: array(0:)
664-
integer(int_size), intent(inout) :: index(0:)
665-
666-
character(len=len(array)) :: tmp
667-
integer(int_size) :: i, tmp_index
668-
669-
tmp = array(0)
670-
tmp_index = index(0)
671-
find_hole: do i=1, size(array, kind=int_size)-1
672-
if ( array(i) >= tmp ) exit find_hole
673-
array(i-1) = array(i)
674-
index(i-1) = index(i)
675-
end do find_hole
676-
array(i-1) = tmp
677-
index(i-1) = tmp_index
678-
679-
end subroutine insert_head
680-
681-
682-
subroutine merge_sort( array, index, buf, ibuf )
683-
! The Rust merge sort borrows some (but not all) of the ideas from TimSort,
684-
! which is described in detail at
685-
! (http://svn.python.org/projects/python/trunk/Objects/listsort.txt).
686-
!
687-
! The algorithm identifies strictly descending and non-descending
688-
! subsequences, which are called natural runs. Where these runs are less
689-
! than a minimum run size they are padded by adding additional samples
690-
! using an insertion sort. The merge process is driven by a stack of
691-
! pending unmerged runs. Each newly found run is pushed onto the stack,
692-
! and then pairs of adjacentd runs are merged until these two invariants
693-
! are satisfied:
694-
!
695-
! 1. for every `i` in `1..size(runs)-1`: `runs(i - 1)%len > runs(i)%len`
696-
! 2. for every `i` in `2..size(runs)-1`: `runs(i - 2)%len >
697-
! runs(i - 1)%len + runs(i)%len`
698-
!
699-
! The invariants ensure that the total running time is `O(n log n)`
700-
! worst-case. Consistency of the indices in `index` with the elements of
701-
! `array` are maintained.
702-
703-
character(len=*), intent(inout) :: array(0:)
704-
integer(int_size), intent(inout) :: index(0:)
705-
character(len=len(array)), intent(inout) :: buf(0:)
706-
integer(int_size), intent(inout) :: ibuf(0:)
707-
708-
integer(int_size) :: array_size, finish, min_run, r, r_count, &
709-
start
710-
type(run_type) :: runs(0:max_merge_stack-1), left, right
711-
712-
array_size = size(array, kind=int_size)
713-
714-
! Very short runs are extended using insertion sort to span at least this
715-
! many elements. Slices of up to this length are sorted using insertion sort.
716-
min_run = calc_min_run( array_size )
717-
718-
if ( array_size <= min_run ) then
719-
if ( array_size >= 2 ) call insertion_sort( array, index )
720-
return
721-
end if
722-
723-
724-
! Following Rust sort, natural runs in `array` are identified by traversing
725-
! it backwards. By traversing it backward, merges more often go in the
726-
! opposite direction (forwards). According to developers of Rust sort,
727-
! merging forwards is slightly faster than merging backwards. Therefore
728-
! identifying runs by traversing backwards should improve performance.
729-
r_count = 0
730-
finish = array_size - 1
731-
do while ( finish >= 0 )
732-
! Find the next natural run, and reverse it if it's strictly descending.
733-
start = finish
734-
if ( start > 0 ) then
735-
start = start - 1
736-
if ( array(start+1) < array(start) ) then
737-
Descending: do while ( start > 0 )
738-
if ( array(start) >= array(start-1) ) &
739-
exit Descending
740-
start = start - 1
741-
end do Descending
742-
call reverse_segment( array(start:finish), &
743-
index(start:finish) )
744-
else
745-
Ascending: do while( start > 0 )
746-
if ( array(start) < array(start-1) ) exit Ascending
747-
start = start - 1
748-
end do Ascending
749-
end if
750-
end if
751-
752-
! If the run is too short insert some more elements using an insertion sort.
753-
Insert: do while ( start > 0 )
754-
if ( finish - start >= min_run - 1 ) exit Insert
755-
start = start - 1
756-
call insert_head( array(start:finish), index(start:finish) )
757-
end do Insert
758-
if ( start == 0 .and. finish == array_size - 1 ) return
759-
760-
runs(r_count) = run_type( base = start, &
761-
len = finish - start + 1 )
762-
finish = start-1
763-
r_count = r_count + 1
764-
765-
! Determine whether pairs of adjacent runs need to be merged to satisfy
766-
! the invariants, and, if so, merge them.
767-
Merge_loop: do
768-
r = collapse( runs(0:r_count - 1) )
769-
if ( r < 0 .or. r_count <= 1 ) exit Merge_loop
770-
left = runs( r + 1 )
771-
right = runs( r )
772-
call merge( array( left % base: &
773-
right % base + right % len - 1 ), &
774-
left % len, buf, &
775-
index( left % base: &
776-
right % base + right % len - 1 ), ibuf )
777-
778-
runs(r) = run_type( base = left % base, &
779-
len = left % len + right % len )
780-
if ( r == r_count - 3 ) runs(r+1) = runs(r+2)
781-
r_count = r_count - 1
782-
783-
end do Merge_loop
784-
end do
785-
if ( r_count /= 1 ) &
786-
error stop "MERGE_SORT completed without RUN COUNT == 1."
787-
788-
end subroutine merge_sort
789-
790-
791-
pure subroutine merge( array, mid, buf, index, ibuf )
792-
! Merges the two non-decreasing runs `ARRAY(0:MID-1)` and `ARRAY(MID:)`
793-
! using `BUF` as temporary storage, and stores the merged runs into
794-
! `ARRAY(0:)`. `MID` must be > 0, and < `SIZE(ARRAY)-1`. Buffer `BUF`
795-
! must be long enough to hold the shorter of the two runs.
796-
character(len=*), intent(inout) :: array(0:)
797-
integer(int_size), intent(in) :: mid
798-
character(len=len(array)), intent(inout) :: buf(0:)
799-
integer(int_size), intent(inout) :: index(0:)
800-
integer(int_size), intent(inout) :: ibuf(0:)
801-
802-
integer(int_size) :: array_len, i, j, k
803-
804-
array_len = size(array, kind=int_size)
805-
806-
! Merge first copies the shorter run into `buf`. Then, depending on which
807-
! run was shorter, it traces the copied run and the longer run forwards
808-
! (or backwards), comparing their next unprocessed elements and then
809-
! copying the lesser (or greater) one into `array`.
810-
811-
if ( mid <= array_len - mid ) then ! The left run is shorter.
812-
buf(0:mid-1) = array(0:mid-1)
813-
ibuf(0:mid-1) = index(0:mid-1)
814-
i = 0
815-
j = mid
816-
merge_lower: do k = 0, array_len-1
817-
if ( buf(i) <= array(j) ) then
818-
array(k) = buf(i)
819-
index(k) = ibuf(i)
820-
i = i + 1
821-
if ( i >= mid ) exit merge_lower
822-
else
823-
array(k) = array(j)
824-
index(k) = index(j)
825-
j = j + 1
826-
if ( j >= array_len ) then
827-
array(k+1:) = buf(i:mid-1)
828-
index(k+1:) = ibuf(i:mid-1)
829-
exit merge_lower
830-
end if
831-
end if
832-
end do merge_lower
833-
else ! The right run is shorter
834-
buf(0:array_len-mid-1) = array(mid:array_len-1)
835-
ibuf(0:array_len-mid-1) = index(mid:array_len-1)
836-
i = mid - 1
837-
j = array_len - mid -1
838-
merge_upper: do k = array_len-1, 0, -1
839-
if ( buf(j) >= array(i) ) then
840-
array(k) = buf(j)
841-
index(k) = ibuf(j)
842-
j = j - 1
843-
if ( j < 0 ) exit merge_upper
844-
else
845-
array(k) = array(i)
846-
index(k) = index(i)
847-
i = i - 1
848-
if ( i < 0 ) then
849-
array(0:k-1) = buf(0:j)
850-
index(0:k-1) = ibuf(0:j)
851-
exit merge_upper
852-
end if
853-
end if
854-
end do merge_upper
855-
end if
856-
end subroutine merge
857-
858-
859-
pure subroutine reverse_segment( array, index )
860-
! Reverse a segment of an array in place
861-
character(len=*), intent(inout) :: array(0:)
862-
integer(int_size), intent(inout) :: index(0:)
863-
864-
integer(int_size) :: itemp, lo, hi
865-
character(len=len(array)) :: temp
866-
867-
lo = 0
868-
hi = size( array, kind=int_size ) - 1
869-
do while( lo < hi )
870-
temp = array(lo)
871-
array(lo) = array(hi)
872-
array(hi) = temp
873-
itemp = index(lo)
874-
index(lo) = index(hi)
875-
index(hi) = itemp
876-
lo = lo + 1
877-
hi = hi - 1
878-
end do
879-
880-
end subroutine reverse_segment
881-
882-
end subroutine char_sort_index
883-
884490
end submodule stdlib_sorting_sort_index

0 commit comments

Comments
 (0)
Please sign in to comment.