This file is subjecâ€ both to the Fortran Standard Library license, and to additional licensing requirements as it contains translations of other software.

The Fortran Standard Library, including this file, is distributed under the MIT license that should be included with the library's distribution.

Copyright (c) 2021 Fortran stdlib developers

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sellcopies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

The generic subroutine, `SORT_INDEX`

, is substantially a translation to
Fortran 2008 of the `"Rust" sort`

sorting routines in
`slice.rs`

The `rust sort`

implementation is distributed with the header:

Copyright 2012-2015 The Rust Project Developers. See the COPYRIGHT file at the top-level directory of this distribution and at http://rust-lang.org/COPYRIGHT.

Licensed under the Apache License, Version 2.0

so the license for the original`slice.rs`

code is compatible with the use
of modified versions of the code in the Fortran Standard Library under
the MIT license.

#:include "common.fypp" #:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_TYPES, INT_KINDS)) #:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_TYPES, REAL_KINDS)) #:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_TYPES, STRING_KINDS)) #:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"], ["character(len=:)"], ["character(len=len(array))"], ["char"])) #! For better code reuse in fypp, make lists that contain the input types, #! with each having output types and a separate name prefix for subroutines #! This approach allows us to have the same code for all input types. #:set IRSC_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME !! Licensing: !! !! This file is subjecâ€ both to the Fortran Standard Library license, and !! to additional licensing requirements as it contains translations of !! other software. !! !! The Fortran Standard Library, including this file, is distributed under !! the MIT license that should be included with the library's distribution. !! !! Copyright (c) 2021 Fortran stdlib developers !! !! Permission is hereby granted, free of charge, to any person obtaining a !! copy of this software and associated documentation files (the !! "Software"), to deal in the Software without restriction, including !! without limitation the rights to use, copy, modify, merge, publish, !! distribute, sublicense, and/or sellcopies of the Software, and to permit !! persons to whom the Software is furnished to do so, subject to the !! following conditions: !! !! The above copyright notice and this permission notice shall be included !! in all copies or substantial portions of the Software. !! !! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS !! OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF !! MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. !! IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY !! CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, !! TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE !! SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. !! !! The generic subroutine, `SORT_INDEX`, is substantially a translation to !! Fortran 2008 of the `"Rust" sort` sorting routines in !! [`slice.rs`](https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs) !! The `rust sort` implementation is distributed with the header: !! !! Copyright 2012-2015 The Rust Project Developers. See the COPYRIGHT !! file at the top-level directory of this distribution and at !! http://rust-lang.org/COPYRIGHT. !! !! Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or !! http://www.apache.org/licenses/LICENSE-2.0> or the MIT license !! <LICENSE-MIT or http://opensource.org/licenses/MIT>, at your !! option. This file may not be copied, modified, or distributed !! except according to those terms. !! !! so the license for the original`slice.rs` code is compatible with the use !! of modified versions of the code in the Fortran Standard Library under !! the MIT license. submodule(stdlib_sorting) stdlib_sorting_sort_index implicit none contains #:for t1, t2, t3, name1 in IRSC_TYPES_ALT_NAME module subroutine ${name1}$_sort_index( array, index, work, iwork, reverse ) ! A modification of `${name1}$_ord_sort` to return an array of indices that ! would perform a stable sort of the `ARRAY` as input, and also sort `ARRAY` ! as desired. The indices by default ! correspond to a non-decreasing sort, but if the optional argument ! `REVERSE` is present with a value of `.TRUE.` the indices correspond to ! a non-increasing sort. The logic of the determination of indexing largely ! follows the `"Rust" sort` found in `slice.rs`: ! https://github.com/rust-lang/rust/blob/90eb44a5897c39e3dff9c7e48e3973671dcd9496/src/liballoc/slice.rs#L2159 ! The Rust version is a simplification of the Timsort algorithm described ! in https://svn.python.org/projects/python/trunk/Objects/listsort.txt, as ! it drops both the use of 'galloping' to identify bounds of regions to be ! sorted and the estimation of the optimal `run size`. However it remains ! a hybrid sorting algorithm combining an iterative Merge sort controlled ! by a stack of `RUNS` identified by regions of uniformly decreasing or ! non-decreasing sequences that may be expanded to a minimum run size, with ! an insertion sort. ! ! Note the Fortran implementation simplifies the logic as it only has to ! deal with Fortran arrays of intrinsic types and not the full generality ! of Rust's arrays and lists for arbitrary types. It also adds the ! estimation of the optimal `run size` as suggested in Tim Peters' ! original listsort.txt, and the optional `work` and `iwork` arraya to be ! used as scratch memory. ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(out) :: index(0:) ${t3}$, intent(out), optional :: work(0:) integer(int_size), intent(out), optional :: iwork(0:) logical, intent(in), optional :: reverse integer(int_size) :: array_size, i, stat ${t2}$, allocatable :: buf(:) integer(int_size), allocatable :: ibuf(:) array_size = size(array, kind=int_size) do i = 0, array_size-1 index(i) = i+1 end do if ( optval(reverse, .false.) ) then call reverse_segment( array, index ) end if ! If necessary allocate buffers to serve as scratch memory. if ( present(work) ) then if ( size(work, kind=int_size) < array_size/2 ) then error stop "work array is too small." end if if ( present(iwork) ) then if ( size(iwork, kind=int_size) < array_size/2 ) then error stop "iwork array is too small." endif call merge_sort( array, index, work, iwork ) else allocate( ibuf(0:array_size/2-1), stat=stat ) if ( stat /= 0 ) error stop "Allocation of index buffer failed." call merge_sort( array, index, work, ibuf ) end if else #:if t1[0:4] == "char" allocate( ${t3}$ :: buf(0:array_size/2-1), & stat=stat ) #:else allocate( buf(0:array_size/2-1), stat=stat ) #:endif if ( stat /= 0 ) error stop "Allocation of array buffer failed." if ( present(iwork) ) then if ( size(iwork, kind=int_size) < array_size/2 ) then error stop "iwork array is too small." endif call merge_sort( array, index, buf, iwork ) else allocate( ibuf(0:array_size/2-1), stat=stat ) if ( stat /= 0 ) error stop "Allocation of index buffer failed." call merge_sort( array, index, buf, ibuf ) end if end if if ( optval(reverse, .false.) ) then call reverse_segment( array, index ) end if contains pure function calc_min_run( n ) result(min_run) !! Returns the minimum length of a run from 32-63 so that N/MIN_RUN is !! less than or equal to a power of two. See !! https://svn.python.org/projects/python/trunk/Objects/listsort.txt integer(int_size) :: min_run integer(int_size), intent(in) :: n integer(int_size) :: num, r num = n r = 0_int_size do while( num >= 64 ) r = ior( r, iand(num, 1_int_size) ) num = ishft(num, -1_int_size) end do min_run = num + r end function calc_min_run pure subroutine insertion_sort( array, index ) ! Sorts `ARRAY` using an insertion sort, while maintaining consistency in ! location of the indices in `INDEX` to the elements of `ARRAY`. ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(inout) :: index(0:) integer(int_size) :: i, j, key_index ${t3}$ :: key do j=1, size(array, kind=int_size)-1 key = array(j) key_index = index(j) i = j - 1 do while( i >= 0 ) if ( array(i) <= key ) exit array(i+1) = array(i) index(i+1) = index(i) i = i - 1 end do array(i+1) = key index(i+1) = key_index end do end subroutine insertion_sort pure function collapse( runs ) result ( r ) ! Examine the stack of runs waiting to be merged, identifying adjacent runs ! to be merged until the stack invariants are restablished: ! ! 1. len(-3) > len(-2) + len(-1) ! 2. len(-2) > len(-1) integer(int_size) :: r type(run_type), intent(in), target :: runs(0:) integer(int_size) :: n logical :: test n = size(runs, kind=int_size) test = .false. if (n >= 2) then if ( runs( n-1 ) % base == 0 .or. & runs( n-2 ) % len <= runs(n-1) % len ) then test = .true. else if ( n >= 3 ) then ! X exists if ( runs(n-3) % len <= & runs(n-2) % len + runs(n-1) % len ) then test = .true. ! |X| <= |Y| + |Z| => will need to merge due to rho1 or rho2 else if( n >= 4 ) then if ( runs(n-4) % len <= & runs(n-3) % len + runs(n-2) % len ) then test = .true. ! |W| <= |X| + |Y| => will need to merge due to rho1 or rho3 end if end if end if end if if ( test ) then ! By default merge Y & Z, rho2 or rho3 if ( n >= 3 ) then if ( runs(n-3) % len < runs(n-1) % len ) then r = n - 3 ! |X| < |Z| => merge X & Y, rho1 return end if end if r = n - 2 ! |Y| <= |Z| => merge Y & Z, rho4 return else r = -1 end if end function collapse pure subroutine insert_head( array, index ) ! Inserts `array(0)` into the pre-sorted sequence `array(1:)` so that the ! whole `array(0:)` becomes sorted, copying the first element into ! a temporary variable, iterating until the right place for it is found. ! copying every traversed element into the slot preceding it, and finally, ! copying data from the temporary variable into the resulting hole. ! Consistency of the indices in `index` with the elements of `array` ! are maintained. ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(inout) :: index(0:) ${t3}$ :: tmp integer(int_size) :: i, tmp_index tmp = array(0) tmp_index = index(0) find_hole: do i=1, size(array, kind=int_size)-1 if ( array(i) >= tmp ) exit find_hole array(i-1) = array(i) index(i-1) = index(i) end do find_hole array(i-1) = tmp index(i-1) = tmp_index end subroutine insert_head subroutine merge_sort( array, index, buf, ibuf ) ! The Rust merge sort borrows some (but not all) of the ideas from TimSort, ! which is described in detail at ! (http://svn.python.org/projects/python/trunk/Objects/listsort.txt). ! ! The algorithm identifies strictly descending and non-descending ! subsequences, which are called natural runs. Where these runs are less ! than a minimum run size they are padded by adding additional samples ! using an insertion sort. The merge process is driven by a stack of ! pending unmerged runs. Each newly found run is pushed onto the stack, ! and then pairs of adjacentd runs are merged until these two invariants ! are satisfied: ! ! 1. for every `i` in `1..size(runs)-1`: `runs(i - 1)%len > runs(i)%len` ! 2. for every `i` in `2..size(runs)-1`: `runs(i - 2)%len > ! runs(i - 1)%len + runs(i)%len` ! ! The invariants ensure that the total running time is `O(n log n)` ! worst-case. Consistency of the indices in `index` with the elements of ! `array` are maintained. ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(inout) :: index(0:) ${t3}$, intent(inout) :: buf(0:) integer(int_size), intent(inout) :: ibuf(0:) integer(int_size) :: array_size, finish, min_run, r, r_count, & start type(run_type) :: runs(0:max_merge_stack-1), left, right array_size = size(array, kind=int_size) ! Very short runs are extended using insertion sort to span at least this ! many elements. Slices of up to this length are sorted using insertion sort. min_run = calc_min_run( array_size ) if ( array_size <= min_run ) then if ( array_size >= 2 ) call insertion_sort( array, index ) return end if ! Following Rust sort, natural runs in `array` are identified by traversing ! it backwards. By traversing it backward, merges more often go in the ! opposite direction (forwards). According to developers of Rust sort, ! merging forwards is slightly faster than merging backwards. Therefore ! identifying runs by traversing backwards should improve performance. r_count = 0 finish = array_size - 1 do while ( finish >= 0 ) ! Find the next natural run, and reverse it if it's strictly descending. start = finish if ( start > 0 ) then start = start - 1 if ( array(start+1) < array(start) ) then Descending: do while ( start > 0 ) if ( array(start) >= array(start-1) ) & exit Descending start = start - 1 end do Descending call reverse_segment( array(start:finish), & index(start:finish) ) else Ascending: do while( start > 0 ) if ( array(start) < array(start-1) ) exit Ascending start = start - 1 end do Ascending end if end if ! If the run is too short insert some more elements using an insertion sort. Insert: do while ( start > 0 ) if ( finish - start >= min_run - 1 ) exit Insert start = start - 1 call insert_head( array(start:finish), index(start:finish) ) end do Insert if ( start == 0 .and. finish == array_size - 1 ) return runs(r_count) = run_type( base = start, & len = finish - start + 1 ) finish = start-1 r_count = r_count + 1 ! Determine whether pairs of adjacent runs need to be merged to satisfy ! the invariants, and, if so, merge them. Merge_loop: do r = collapse( runs(0:r_count - 1) ) if ( r < 0 .or. r_count <= 1 ) exit Merge_loop left = runs( r + 1 ) right = runs( r ) call merge( array( left % base: & right % base + right % len - 1 ), & left % len, buf, & index( left % base: & right % base + right % len - 1 ), ibuf ) runs(r) = run_type( base = left % base, & len = left % len + right % len ) if ( r == r_count - 3 ) runs(r+1) = runs(r+2) r_count = r_count - 1 end do Merge_loop end do if ( r_count /= 1 ) & error stop "MERGE_SORT completed without RUN COUNT == 1." end subroutine merge_sort pure subroutine merge( array, mid, buf, index, ibuf ) ! Merges the two non-decreasing runs `ARRAY(0:MID-1)` and `ARRAY(MID:)` ! using `BUF` as temporary storage, and stores the merged runs into ! `ARRAY(0:)`. `MID` must be > 0, and < `SIZE(ARRAY)-1`. Buffer `BUF` ! must be long enough to hold the shorter of the two runs. ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(in) :: mid ${t3}$, intent(inout) :: buf(0:) integer(int_size), intent(inout) :: index(0:) integer(int_size), intent(inout) :: ibuf(0:) integer(int_size) :: array_len, i, j, k array_len = size(array, kind=int_size) ! Merge first copies the shorter run into `buf`. Then, depending on which ! run was shorter, it traces the copied run and the longer run forwards ! (or backwards), comparing their next unprocessed elements and then ! copying the lesser (or greater) one into `array`. if ( mid <= array_len - mid ) then ! The left run is shorter. buf(0:mid-1) = array(0:mid-1) ibuf(0:mid-1) = index(0:mid-1) i = 0 j = mid merge_lower: do k = 0, array_len-1 if ( buf(i) <= array(j) ) then array(k) = buf(i) index(k) = ibuf(i) i = i + 1 if ( i >= mid ) exit merge_lower else array(k) = array(j) index(k) = index(j) j = j + 1 if ( j >= array_len ) then array(k+1:) = buf(i:mid-1) index(k+1:) = ibuf(i:mid-1) exit merge_lower end if end if end do merge_lower else ! The right run is shorter buf(0:array_len-mid-1) = array(mid:array_len-1) ibuf(0:array_len-mid-1) = index(mid:array_len-1) i = mid - 1 j = array_len - mid -1 merge_upper: do k = array_len-1, 0, -1 if ( buf(j) >= array(i) ) then array(k) = buf(j) index(k) = ibuf(j) j = j - 1 if ( j < 0 ) exit merge_upper else array(k) = array(i) index(k) = index(i) i = i - 1 if ( i < 0 ) then array(0:k-1) = buf(0:j) index(0:k-1) = ibuf(0:j) exit merge_upper end if end if end do merge_upper end if end subroutine merge pure subroutine reverse_segment( array, index ) ! Reverse a segment of an array in place ${t1}$, intent(inout) :: array(0:) integer(int_size), intent(inout) :: index(0:) integer(int_size) :: itemp, lo, hi ${t3}$ :: temp lo = 0 hi = size( array, kind=int_size ) - 1 do while( lo < hi ) temp = array(lo) array(lo) = array(hi) array(hi) = temp itemp = index(lo) index(lo) = index(hi) index(hi) = itemp lo = lo + 1 hi = hi - 1 end do end subroutine reverse_segment end subroutine ${name1}$_sort_index #:endfor end submodule stdlib_sorting_sort_index