pure subroutine sp_increase_sort(array)
! `sp_increase_sort( array )` sorts the input `ARRAY` of type `real(sp)`
! using a hybrid sort based on the `introsort` of David Musser. As with
! `introsort`, `sp_increase_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.
real(sp), 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`.
real(sp), 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.
real(sp), intent(inout) :: array(0:)
integer(int_index), intent(out) :: index
real(sp) :: 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 > v) .neqv. (u > 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 < u) .neqv. (v < 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) <= 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.
real(sp), intent(inout) :: array(0:)
integer(int_index) :: i, j
real(sp) :: 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) <= 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
real(sp), intent(inout) :: array(0:)
integer(int_index) :: i, heap_size
real(sp) :: 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
real(sp), intent(inout) :: array(0:)
integer(int_index), intent(in) :: i, heap_size
integer(int_index) :: l, r, largest
real(sp) :: y
largest = i
l = 2_int_index*i + 1_int_index
r = l + 1_int_index
if (l < heap_size) then
if (array(l) > array(largest)) largest = l
end if
if (r < heap_size) then
if (array(r) > 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 sp_increase_sort