stdlib_sorting_sort.fypp Source File

#:include "common.fypp"
#:set INT_TYPES_ALT_NAME = list(zip(INT_TYPES, INT_TYPES, INT_KINDS))
#:set REAL_TYPES_ALT_NAME = list(zip(REAL_TYPES, REAL_TYPES, REAL_KINDS))
#:set STRING_TYPES_ALT_NAME = list(zip(STRING_TYPES, STRING_TYPES, STRING_KINDS))
#:set CHAR_TYPES_ALT_NAME = list(zip(["character(len=*)"],  ["character(len=len(array))"], ["char"]))
#:set BITSET_TYPES_ALT_NAME = list(zip(BITSET_TYPES, BITSET_TYPES, BITSET_KINDS))

#! 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 IRSCB_TYPES_ALT_NAME = INT_TYPES_ALT_NAME + REAL_TYPES_ALT_NAME + STRING_TYPES_ALT_NAME + CHAR_TYPES_ALT_NAME &
    & + BITSET_TYPES_ALT_NAME

#: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 t1, t2, name1 in IRSCB_TYPES_ALT_NAME
    pure module subroutine ${name1}$_sort( array, reverse )
            ${t1}$, intent(inout) :: array(0:)
            logical, intent(in), optional            :: reverse

            if(optval(reverse, .false.))then
             call ${name1}$_decrease_sort(array)
            else
             call ${name1}$_increase_sort(array)
            endif
    end subroutine ${name1}$_sort
#:endfor

#:for sname, signt, signoppt in SIGN_NAME_TYPE
#:for t1, t2, name1 in IRSCB_TYPES_ALT_NAME

    pure subroutine ${name1}$_${sname}$_sort( array )
! `${name1}$_${sname}$_sort( array )` sorts the input `ARRAY` of type `${t1}$`
! using a hybrid sort based on the `introsort` of David Musser. As with
! `introsort`, `${name1}$_${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_index),  &
                                                 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_index), parameter :: insert_size = 16_int_index
            integer(int_index)            :: index

            if ( size(array, kind=int_index) <= 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_index), intent(out) :: index

            ${t2}$ :: u, v, w, x, y
            integer(int_index) :: i, j

! Determine median of three and exchange it with the end.
            u = array( 0 )
            v = array( size(array, kind=int_index)/2-1 )
            w = array( size(array, kind=int_index)-1 )
            if ( (u ${signt}$ v) .neqv. (u ${signt}$ w) ) then
                x = u
                y = array(0)
                array(0) = array( size( array, kind=int_index ) - 1 )
                array( size( array, kind=int_index ) - 1 ) = y
            else if ( (v ${signoppt}$ u) .neqv. (v ${signoppt}$ w) ) then
                x = v
                y = array(size( array, kind=int_index )/2-1)
                array( size( array, kind=int_index )/2-1 ) = &
                    array( size( array, kind=int_index )-1 )
                array( size( array, kind=int_index )-1 ) = y
            else
                x = w
            end if
! Partition the array.
            i = -1_int_index
            do j = 0_int_index, size(array, kind=int_index)-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_index)-1)
            array(size(array, kind=int_index)-1) = y
            index = i + 1

        end subroutine partition

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

            integer(int_index) :: i, j
            ${t2}$ :: key

            do j=1_int_index, size(array, kind=int_index)-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_index) :: i, heap_size
            ${t2}$   :: y

            heap_size = size( array, kind=int_index )
! Build the max heap
            do i = (heap_size-2)/2_int_index, 0_int_index, -1_int_index
                call max_heapify( array, i, heap_size )
            end do
            do i = heap_size-1, 1_int_index, -1_int_index
! 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_index, 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_index), intent(in)  :: i, heap_size

            integer(int_index) :: l, r, largest
            ${t2}$   :: y

            largest = i
            l = 2_int_index * i + 1_int_index
            r = l + 1_int_index
            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 ${name1}$_${sname}$_sort

#:endfor
#:endfor

end submodule stdlib_sorting_sort