stdlib_sorting_sort.fypp

#:include "common.fypp"
#:set IRS_KINDS_TYPES = INT_KINDS_TYPES + REAL_KINDS_TYPES + STRING_KINDS_TYPES

#:set SIGN_NAME = ["increase", "decrease"]
#:set SIGN_TYPE = [">", "<"]
#:set SIGN_OPP_TYPE = ["<", ">"]
#:set SIGN_NAME_TYPE = list(zip(SIGN_NAME, SIGN_TYPE, SIGN_OPP_TYPE))

!! 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`, is substantially a
!! translation to Fortran 2008, of the `introsort` of David Musser.
!! David Musser has given permission to include a variant of `introsort`
!! in the Fortran Standard Library under the MIT license provided
!! we cite:
!!
!!   Musser, D.R., “Introspective Sorting and Selection Algorithms,”
!!   Software—Practice and Experience, Vol. 27(8), 983–993 (August 1997).
!!
!! as the official source of the algorithm.

submodule(stdlib_sorting) stdlib_sorting_sort
!! This submodule implements the overloaded sorting subroutine `SORT`
!! that can be used to sort four kinds of `INTEGER` arrays and three kinds
!! of `REAL` arrays. Sorting is in order of increasing value, with the worst
!! case run time performance of `O(N Ln(N))`.
!!
!! `SORT` uses the `INTROSORT` sorting algorithm of David Musser,
!! http://www.cs.rpi.edu/~musser/gp/introsort.ps. `introsort` is a hybrid
!! unstable comparison algorithm combining `quicksort`, `insertion sort`, and
!! `heap sort`. While this algorithm is always O(N Ln(N)) it is relatively
!! fast on randomly ordered data, but inconsistent in performance on partly
!! sorted data, sometimes having `merge sort` performance, sometimes having
!! better than `quicksort` performance.

    implicit none

contains

#:for k1, t1 in IRS_KINDS_TYPES
    pure module subroutine ${k1}$_sort( array, reverse )
            ${t1}$, intent(inout) :: array(0:)
            logical, intent(in), optional            :: reverse

            logical :: reverse_

            reverse_ = .false.
            if(present(reverse)) reverse_ = reverse

            if(reverse_)then
             call ${k1}$_decrease_sort(array)
            else
             call ${k1}$_increase_sort(array)
            endif
    end subroutine ${k1}$_sort
#:endfor

#:for sname, signt, signoppt in SIGN_NAME_TYPE
#:for k1, t1 in IRS_KINDS_TYPES

    pure subroutine ${k1}$_${sname}$_sort( array )
! `${k1}$_${sname}$_sort( array )` sorts the input `ARRAY` of type `${t1}$`
! using a hybrid sort based on the `introsort` of David Musser. As with
! `introsort`, `${k1}$_${sname}$_sort( array )` is an unstable hybrid comparison
! algorithm using `quicksort` for the main body of the sort tree,
! supplemented by `insertion sort` for the outer branches, but if
! `quicksort` is converging too slowly the algorithm resorts
! to `heapsort`. The algorithm is of order O(N Ln(N)) for all inputs.
! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
! behavior is typically small compared to other sorting algorithms.

        ${t1}$, intent(inout) :: array(0:)

        integer(int32) :: depth_limit

        depth_limit = 2 * int( floor( log( real( size( array, kind=int_size),  &
                                                 kind=dp) ) / log(2.0_dp) ), &
                               kind=int32 )
        call introsort(array, depth_limit)

    contains

        pure recursive subroutine introsort( array, depth_limit )
! It devolves to `insertionsort` if the remaining number of elements
! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion
! of the `quicksort` is too slow as estimated from `DEPTH_LIMIT`,
! otherwise sorting is done by a `quicksort`.
            ${t1}$, intent(inout) :: array(0:)
            integer(int32), intent(in)     :: depth_limit

            integer(int_size), parameter :: insert_size = 16_int_size
            integer(int_size)            :: index

            if ( size(array, kind=int_size) <= insert_size ) then
                ! May be best at the end of SORT processing the whole array
                ! See Musser, D.R., “Introspective Sorting and Selection
                ! Algorithms,” Software—Practice and Experience, Vol. 27(8),
                ! 983–993 (August 1997).

                call insertion_sort( array )
            else if ( depth_limit == 0 ) then
                call heap_sort( array )
            else
                call partition( array, index )
                call introsort( array(0:index-1), depth_limit-1 )
                call introsort( array(index+1:), depth_limit-1 )
            end if

        end subroutine introsort


        pure subroutine partition( array, index )
! quicksort partition using median of three.
            ${t1}$, intent(inout) :: array(0:)
            integer(int_size), intent(out) :: index

            ${t1}$ :: u, v, w, x, y
            integer(int_size) :: i, j

! Determine median of three and exchange it with the end.
            u = array( 0 )
            v = array( size(array, kind=int_size)/2-1 )
            w = array( size(array, kind=int_size)-1 )
            if ( (u ${signt}$ v) .neqv. (u ${signt}$ w) ) then
                x = u
                y = array(0)
                array(0) = array( size( array, kind=int_size ) - 1 )
                array( size( array, kind=int_size ) - 1 ) = y
            else if ( (v ${signoppt}$ u) .neqv. (v ${signoppt}$ w) ) then
                x = v
                y = array(size( array, kind=int_size )/2-1)
                array( size( array, kind=int_size )/2-1 ) = &
                    array( size( array, kind=int_size )-1 )
                array( size( array, kind=int_size )-1 ) = y
            else
                x = w
            end if
! Partition the array.
            i = -1_int_size
            do j = 0_int_size, size(array, kind=int_size)-2
                if ( array(j) ${signoppt}$= x ) then
                    i = i + 1
                    y = array(i)
                    array(i) = array(j)
                    array(j) = y
                end if
            end do
            y = array(i+1)
            array(i+1) = array(size(array, kind=int_size)-1)
            array(size(array, kind=int_size)-1) = y
            index = i + 1

        end subroutine partition

        pure subroutine insertion_sort( array )
! Bog standard insertion sort.
            ${t1}$, intent(inout) :: array(0:)

            integer(int_size) :: i, j
            ${t1}$ :: key

            do j=1_int_size, size(array, kind=int_size)-1
                key = array(j)
                i = j - 1
                do while( i >= 0 )
                    if ( array(i) ${signoppt}$= key ) exit
                    array(i+1) = array(i)
                    i = i - 1
                end do
                array(i+1) = key
            end do

        end subroutine insertion_sort

        pure subroutine heap_sort( array )
! A bog standard heap sort
            ${t1}$, intent(inout) :: array(0:)

            integer(int_size) :: i, heap_size
            ${t1}$   :: y

            heap_size = size( array, kind=int_size )
! Build the max heap
            do i = (heap_size-2)/2_int_size, 0_int_size, -1_int_size
                call max_heapify( array, i, heap_size )
            end do
            do i = heap_size-1, 1_int_size, -1_int_size
! Swap the first element with the current final element
                y = array(0)
                array(0) = array(i)
                array(i) = y
! Sift down using max_heapify
                call max_heapify( array, 0_int_size, i )
            end do

        end subroutine heap_sort

        pure recursive subroutine max_heapify( array, i, heap_size )
! Transform the array into a max heap
            ${t1}$, intent(inout) :: array(0:)
            integer(int_size), intent(in)  :: i, heap_size

            integer(int_size) :: l, r, largest
            ${t1}$   :: y

            largest = i
            l = 2_int_size * i + 1_int_size
            r = l + 1_int_size
            if ( l < heap_size ) then
                if ( array(l) ${signt}$ array(largest) ) largest = l
            end if
            if ( r < heap_size ) then
                if ( array(r) ${signt}$ array(largest) ) largest = r
            end if
            if ( largest /= i ) then
                y = array(i)
                array(i) = array(largest)
                array(largest) = y
                call max_heapify( array, largest, heap_size )
            end if

        end subroutine max_heapify

    end subroutine ${k1}$_${sname}$_sort

#:endfor
#:endfor


    pure module subroutine char_sort( array, reverse )
        character(len=*), intent(inout) :: array(0:)
        logical, intent(in), optional   :: reverse

        logical :: reverse_

        reverse_ = .false.
        if(present(reverse)) reverse_ = reverse

        if(reverse_)then
            call char_decrease_sort(array)
        else
            call char_increase_sort(array)
        endif
    end subroutine char_sort



#:for sname, signt, signoppt in SIGN_NAME_TYPE
    pure subroutine char_${sname}$_sort( array )
! `char_${sname}$_sort( array )` sorts the input `ARRAY` of type `CHARACTER(*)`
! using a hybrid sort based on the `introsort` of David Musser. As with
! `introsort`, `char_${sname}$_sort( array )` is an unstable hybrid comparison
! algorithm using `quicksort` for the main body of the sort tree,
! supplemented by `insertion sort` for the outer branches, but if
! `quicksort` is converging too slowly the algorithm resorts
! to `heapsort`. The algorithm is of order O(N Ln(N)) for all inputs.
! Because it relies on `quicksort`, the coefficient of the O(N Ln(N))
! behavior is typically small compared to other sorting algorithms.

        character(len=*), intent(inout) :: array(0:)

        integer(int32) :: depth_limit

        depth_limit = 2 * int( floor( log( real( size( array, kind=int_size ),  &
                                                 kind=dp) ) / log(2.0_dp) ), &
                               kind=int32 )
        call introsort(array, depth_limit)

    contains

        pure recursive subroutine introsort( array, depth_limit )
! It devolves to `insertionsort` if the remaining number of elements
! is fewer than or equal to `INSERT_SIZE`, `heapsort` if the completion
! of the `quicksort` is too slow as estimated from `DEPTH_LIMIT`,
! otherwise sorting is done by a `quicksort`.
            character(len=*), intent(inout) :: array(0:)
            integer(int32), intent(in)      :: depth_limit

            integer(int_size), parameter :: insert_size = 16_int_size
            integer(int_size)            :: index

            if ( size(array, kind=int_size) <= insert_size ) then
                ! May be best at the end of SORT processing the whole array
                ! See Musser, D.R., “Introspective Sorting and Selection
                ! Algorithms,” Software—Practice and Experience, Vol. 27(8),
                ! 983–993 (August 1997).

                call insertion_sort( array )
            else if ( depth_limit == 0 ) then
                call heap_sort( array )
            else
                call partition( array, index )
                call introsort( array(0:index-1), depth_limit-1 )
                call introsort( array(index+1:), depth_limit-1 )
            end if

        end subroutine introsort


        pure subroutine partition( array, index )
! quicksort partition using median of three.
            character(len=*), intent(inout) :: array(0:)
            integer(int_size), intent(out)  :: index

            integer(int_size)         :: i, j
            character(len=len(array)) :: u, v, w, x, y

! Determine median of three and exchange it with the end.
            u = array( 0 )
            v = array( size(array, kind=int_size)/2-1 )
            w = array( size(array, kind=int_size)-1 )
            if ( (u ${signt}$ v) .neqv. (u ${signt}$ w) ) then
                x = u
                y = array(0)
                array(0) = array( size( array, kind=int_size ) - 1 )
                array( size( array, kind=int_size ) - 1 ) = y
            else if ( (v ${signoppt}$ u) .neqv. (v ${signoppt}$ w) ) then
                x = v
                y = array(size( array, kind=int_size )/2-1)
                array( size( array, kind=int_size )/2-1 ) = &
                    array( size( array, kind=int_size )-1 )
                array( size( array, kind=int_size )-1 ) = y
            else
                x = w
            end if
! Partition the array.
            i = -1_int_size
            do j = 0_int_size, size(array, kind=int_size)-2
                if ( array(j) ${signoppt}$= x ) then
                    i = i + 1
                    y = array(i)
                    array(i) = array(j)
                    array(j) = y
                end if
            end do
            y = array(i+1)
            array(i+1) = array(size(array, kind=int_size)-1)
            array(size(array, kind=int_size)-1) = y
            index = i + 1

        end subroutine partition

        pure subroutine insertion_sort( array )
! Bog standard insertion sort.
            character(len=*), intent(inout) :: array(0:)

            integer(int_size)         :: i, j
            character(len=len(array)) :: key

            do j=1_int_size, size(array, kind=int_size)-1
                key = array(j)
                i = j - 1
                do while( i >= 0 )
                    if ( array(i) ${signoppt}$= key ) exit
                    array(i+1) = array(i)
                    i = i - 1
                end do
                array(i+1) = key
            end do

        end subroutine insertion_sort

        pure subroutine heap_sort( array )
! A bog standard heap sort
            character(len=*), intent(inout) :: array(0:)

            integer(int_size)         :: i, heap_size
            character(len=len(array)) :: y

            heap_size = size( array, kind=int_size )
! Build the max heap
            do i = (heap_size-2)/2_int_size, 0_int_size, -1_int_size
                call max_heapify( array, i, heap_size )
            end do
            do i = heap_size-1, 1_int_size, -1_int_size
! Swap the first element with the current final element
                y = array(0)
                array(0) = array(i)
                array(i) = y
! Sift down using max_heapify
                call max_heapify( array, 0_int_size, i )
            end do

        end subroutine heap_sort

        pure recursive subroutine max_heapify( array, i, heap_size )
! Transform the array into a max heap
            character(len=*), intent(inout) :: array(0:)
            integer(int_size), intent(in)   :: i, heap_size

            integer(int_size)         :: l, r, largest
            character(len=len(array)) :: y

            largest = i
            l = 2_int_size * i + 1_int_size
            r = l + 1_int_size
            if ( l < heap_size ) then
                if ( array(l) ${signt}$ array(largest) ) largest = l
            end if
            if ( r < heap_size ) then
                if ( array(r) ${signt}$ array(largest) ) largest = r
            end if
            if ( largest /= i ) then
                y = array(i)
                array(i) = array(largest)
                array(largest) = y
                call max_heapify( array, largest, heap_size )
            end if

        end subroutine max_heapify

    end subroutine char_${sname}$_sort
#:endfor

end submodule stdlib_sorting_sort