sp_radix_sort Module Subroutine

   module subroutine sp_radix_sort(array, work, reverse, err)
      use, intrinsic :: iso_c_binding, only: c_loc, c_f_pointer
      real(kind=sp), dimension(:), intent(inout), target :: array
      real(kind=sp), dimension(:), intent(inout), target, optional :: work
      logical, intent(in), optional :: reverse
      type(error_t), intent(inout), optional :: err
      integer(kind=int_index) :: i, N, pos, rev_pos
      integer(kind=int32), dimension(:), pointer :: arri32
      integer(kind=int32), dimension(:), pointer :: buffer
      real(kind=sp) :: 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, &
                                "sp_radix_sort: work array is too small.")
               return
            end if
            error stop "sp_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.
         call c_f_pointer(c_loc(work), buffer, [N])
      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, &
                                "sp_radix_sort: allocation of scratch buffer failed.")
               return
            end if
            error stop "sp_radix_sort: Allocation of buffer failed."
         end if
      end if
      call c_f_pointer(c_loc(array), arri32, [N])
      call radix_sort_u32_helper(N, arri32, buffer)
! After calling `radix_sort_u<width>_helper. The array is sorted as unsigned integers.
! The positive real number is sorted, guaranteed by IEEE-754 standard.
! But the negative real number is sorted in a reversed order, and also in the tail of array.
! Remark that -0.0 is the minimum nagative integer, so using radix sort, -0.0 is naturally lesser than 0.0.
! In IEEE-754 standard, the bit representation of `Inf` is greater than all positive real numbers,
! and the `-Inf` is lesser than all negative real numbers. So the order of `Inf, -Inf` is ok.
! The bit representation of `NaN` may be positive or negative integers in different machine,
! thus if the array contains `NaN`, the result is undefined.
      if (arri32(1) >= 0 .and. arri32(N) < 0) then
         pos = 1
         rev_pos = N
         do while (arri32(rev_pos) < 0)
            buffer(pos) = arri32(rev_pos)
            pos = pos + 1
            rev_pos = rev_pos - 1
         end do
         buffer(pos:N) = arri32(1:rev_pos)
         arri32(:) = 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 sp_radix_sort