Compiler Explorer

Source code

module sorting_f
   use iso_fortran_env, only: int32,int64
   implicit none
   private

integer, parameter, public :: IKIND = int32
   integer, parameter, public :: ISIZE = int64
  
   public :: argsort_f
contains

!> Argsort: Federico's implementation
   subroutine argsort_f(n,list,keys)
      integer(ISIZE), intent(in)    :: n
      integer(ISIZE), intent(inout) :: list(*) ! indices to the sorted array
      integer(IKIND), intent(in)    :: keys(*) ! Array to be sorted

integer :: sorted(n)

sorted = keys(1:n)
      call int_quick_sort_andlist(sorted,list(1:n))

end subroutine argsort_f

pure recursive subroutine int_quick_sort_andlist(list,ilist,down)

integer(IKIND), dimension(:), intent(inout) :: list
      integer(ISIZE), dimension(size(list)), intent(inout) :: ilist
      logical, optional, intent(in) :: down

integer(ISIZE)     :: i, j, n
      integer(ISIZE), parameter :: max_simple_sort_size = 8
      integer(IKIND)     :: chosen
      logical            :: descending

descending = .false.; if (present(down)) descending = down
      n = size(list)

choose_sorting_algorithm: if (n <= max_simple_sort_size) then
         ! Use interchange sort for small lists
         do i = 1, size(list) - 1
            do j = i + 1, size(list)
               if (toBeSwapped(list(i),list(j),.false.)) then
                   call swap(list(i),list(j))
                   call swaplist(ilist(i),ilist(j))
               end if
            end do
         end do

else
         ! Use partition (quick) sort if the list is big
         chosen = list(int(n/2))
         i = 0
         j = n + 1

scan_lists: do
            ! Scan list from left end
            ! until element >= chosen is found
            scan_from_left: do
              i = i + 1
              if (toBeSwapped(list(i),chosen,.true.) .or. i==n) exit scan_from_left
            end do scan_from_left

! Scan list from right end
            ! until element <= chosen is found

scan_from_right: do
               j = j - 1
               if (toBeSwapped(chosen,list(j),.true.) .or. j==1) exit scan_from_right
            end do scan_from_right

swap_element: if (i < j) then
                ! Swap two out of place elements
                call swap(list(i),list(j))
                call swaplist(ilist(i),ilist(j))
            else if (i == j) then
                i = i + 1
                exit
            else
                exit
            endif swap_element

end do scan_lists

if (1 < j) call int_quick_sort_andlist(list(:j),ilist(:j),down)
         if (i < n) call int_quick_sort_andlist(list(i:),ilist(i:),down)

end if choose_sorting_algorithm ! test for small array

contains

elemental logical function toBeSwapped(a,b,orEqual)
            integer(IKIND), intent(in) :: a,b
            logical, intent(in) :: orEqual

toBeSwapped = merge(a<b,a>b,descending)
            if (orEqual .and. a==b) toBeSwapped = .true.

end function toBeSwapped

elemental subroutine swap(a,b)
            integer(IKIND), intent(inout) :: a,b
            integer(IKIND) :: tmp
            tmp = a
            a = b
            b = tmp
         end subroutine swap
         elemental subroutine swaplist(a,b)
            integer(ISIZE), intent(inout) :: a,b
            integer(ISIZE) :: tmp
            tmp = a
            a = b
            b = tmp
         end subroutine swaplist

end subroutine int_quick_sort_andlist

end module sorting_f

module sorting
   use, intrinsic :: iso_c_binding
   implicit none
   private
   public :: argsort
contains

!> Returns the indices that would sort an array
   !>
   !> Arguments:
   !>   n: number of elements
   !>   list: on entry, list of indices 1..(1)..n 
   !>           (for reverse sort, n..(-1)..1)
   !>         on return, the indices such that keys(list) would be sorted
   !>   keys: array to be sorted
   !>
   subroutine argsort(n,list,keys)
      integer, intent(in) :: n
      integer, intent(inout), target :: list(*)
      integer, intent(in) :: keys(*)

interface
         subroutine qsort(ptr,count,size,comp) bind(c,name="qsort")
            import c_ptr, c_int, c_size_t, c_funptr
            implicit none
            type(c_ptr), intent(in), value :: ptr
               !> pointer to the array to sort
            integer(c_size_t), intent(in), value :: count
               !> number of elements in the array
            integer(c_size_t), intent(in), value :: size
               !> size of each element in the array in bytes
            type(c_funptr), intent(in), value :: comp
               !> comparison function which returns a negative integer value 
               !> if the first argument is less than the second, a positive 
               !> integer value if the first argument is greater than the 
               !> second and zero if the arguments are equivalent
         end subroutine
      end interface

if (n < 2) then
         ! list already sorted
         return
      end if

call qsort(c_loc(list),int(n,c_size_t),c_sizeof(list(1)), &
         c_funloc(argsort_compare))

contains
      ! technically, the inputs should be void pointers
      ! but works nevertheless
      integer(c_int) function argsort_compare(a,b) bind(c)
         integer(c_int), intent(in) :: a, b
         argsort_compare = keys(a) - keys(b)
      end function
   end subroutine

end module sorting

program main
   use iso_c_binding, only: c_int
   use sorting
   use sorting_f
   implicit none

integer(ISIZE), parameter :: N = 100000

real, allocatable :: x(:)
   integer(ISIZE), allocatable :: listf(:)
   integer(IKIND), allocatable :: keys(:),list(:)
   integer(ISIZE) :: j
   real :: t0,t1,t2

allocate(list(N),keys(N),x(N))
   call random_number(x)
   keys(:) = nint(N*x)
   list(:) = [(j,j=1,N)]
   listF   = list

call cpu_time(t0); call argsort(int(N,c_int),list,keys)
   print *, keys(1:5),list(1:5)
   call cpu_time(t1); call argsort_f(N,listf,keys)
   print *, keys(1:5),list(1:5)
   call cpu_time(t2)

! Supposed to print [ 3, 5, 6, 7]
   print *, 'C time=',t1-t0
   print *, 'F time=',t2-t1

end program