pure module subroutine int64_radix_sort(array, work, reverse, err)
integer(kind=int64), dimension(:), intent(inout) :: array
integer(kind=int64), dimension(:), intent(inout), target, optional :: work
logical, intent(in), optional :: reverse
type(error_t), intent(inout), optional :: err
integer(kind=int_index) :: i, N, start, middle, end
integer(kind=int64), dimension(:), pointer :: buffer
integer(kind=int64) :: item
logical :: use_internal_buffer
integer :: stat
N = size(array, kind=int_index)
if (present(work)) then
if (size(work, kind=int_index) < N) then
if (present(err)) then
call error_raise(err, ERROR_VALIDATION, &
"int64_radix_sort: work array is too small.")
return
end if
error stop "int64_radix_sort: work array is too small."
end if
end if
! An array of fewer than two elements is already sorted, and reversing one is a
! no-op, so there is nothing left to do once `work` has been validated. This
! return is not an optimisation. The sign rotation further down reads array(1)
! and array(N), which at N == 0 are both outside the array; with bounds checking
! that is a hard error, and without it the two reads pick up whatever happens to
! sit next to the array. If they satisfy `array(1) >= 0 .and. array(N) < 0` the
! rotation enters its binary search with start = 1, end = 0 and middle = 0, and
! every iteration recomputes middle = (1 + 0)/2 = 0 and start = 1, so
! `start == end` is never true and the `do while (.true.)` never exits. That is
! a hang whose occurrence depends on adjacent memory, which is why it appears on
! some platforms and not others. The sp and dp specialisations additionally form
! c_loc on a zero-sized array, which is not legal either.
if (N < 2) return
if (present(work)) then
use_internal_buffer = .false.
buffer => work
else
use_internal_buffer = .true.
! `buffer` is a pointer, so `associated` is what `allocated` is for an
! allocatable; `nullify` first because a failed `allocate` leaves a pointer's
! association status unchanged, i.e. undefined for a fresh local pointer.
nullify (buffer)
allocate (buffer(N), stat=stat)
if (.not. associated(buffer)) then
if (present(err)) then
call error_raise(err, ERROR_ALLOC, &
"int64_radix_sort: allocation of scratch buffer failed.")
return
end if
error stop "int64_radix_sort: Allocation of buffer failed."
end if
end if
call radix_sort_u64_helper(N, array, buffer)
if (array(1) >= 0 .and. array(N) < 0) then
start = 1
end = N
middle = (1 + N)/2
do while (.true.)
if (array(middle) >= 0) then
start = middle + 1
else
end = middle
end if
middle = (start + end)/2
if (start == end) exit
end do
buffer(1:(N - middle + 1)) = array(middle:N)
buffer(N - middle + 2:N) = array(1:middle - 1)
array(:) = buffer(:)
end if
if (pic_optional(reverse, .false.)) then
do i = 1, N/2
item = array(i)
array(i) = array(N - i + 1)
array(N - i + 1) = item
end do
end if
if (use_internal_buffer) then
deallocate (buffer)
end if
end subroutine int64_radix_sort