Module M_sort
use, intrinsic :: ISO_FORTRAN_ENV, only : INT8, INT16, INT32, INT64, REAL32, REAL64, REAL128 
implicit none      
private           
public sort_quick_rx

integer,parameter :: ASCII=kind('A')

interface sort_quick_rx

   module procedure sort_quick_rx_real_real32_int32
   module procedure sort_quick_rx_real_real64_int32
   module procedure sort_quick_rx_integer_int8_int32
   module procedure sort_quick_rx_integer_int16_int32
   module procedure sort_quick_rx_integer_int32_int32
   module procedure sort_quick_rx_integer_int64_int32
   module procedure sort_quick_rx_character_ascii_int32

   module procedure sort_quick_rx_real_real32_int64
   module procedure sort_quick_rx_real_real64_int64
   module procedure sort_quick_rx_integer_int8_int64
   module procedure sort_quick_rx_integer_int16_int64
   module procedure sort_quick_rx_integer_int32_int64
   module procedure sort_quick_rx_integer_int64_int64
   module procedure sort_quick_rx_character_ascii_int64

end interface
contains
$!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$PARCEL SORT_QUICK_RX
subroutine sort_quick_rx_${TYPE}_${KIND}_${IKIND}(data,indx)

$@(#) M_sort::sort_quick_rx_${TYPE}_${KIND}_${IKIND}(3f): indexed hybrid quicksort of a ${TYPE}(kind=${KIND}) array

$IFDEF CHARACTER
${TYPE}(kind=${KIND},len=*),intent(in)   :: data(:)
integer(kind=${IKIND}),intent(out)       :: indx(:)
${TYPE}(kind=${KIND},len=len(data))  :: datap
$ELSE
${TYPE}(kind=${KIND}),intent(in)   :: data(:)
integer(kind=${IKIND}),intent(out) :: indx(:)
${TYPE}(kind=${KIND})  :: datap
$ENDIF

integer(kind=int64)    :: n
integer(kind=int64)    :: lstk(31),rstk(31),istk
integer(kind=int64)    :: l,r,i,j,p,indexp,indext

integer,parameter :: M=9
n=size(data)
if(size(indx).lt.n)then  ! if index is not big enough, only sort part of the data
  write(*,*)'*sort_quick_rx_${TYPE}_${KIND}* ERROR: insufficient space to store index data'
  n=size(indx)
endif

do i=1,n
   indx(i)=i
enddo

if (N.gt.M)then

   istk=0
   l=1
   r=n
   TOP: do

      i=l
      j=r

      p=(l+r)/2
      indexp=indx(p)
      datap=data(indexp)

      if (data(indx(l)) .gt. datap) then
         indx(p)=indx(l)
         indx(l)=indexp
         indexp=indx(p)
         datap=data(indexp)
      endif

      if (datap .gt. data(indx(r))) then

         if (data(indx(l)) .gt. data(indx(r))) then
            indx(p)=indx(l)
            indx(l)=indx(r)
         else
            indx(p)=indx(r)
         endif

         indx(r)=indexp
         indexp=indx(p)
         datap=data(indexp)
      endif

      Q3: do
         I=I+1
         if (data(indx(i)).lt.datap)then
            cycle Q3
         endif
         Q4: do
            j=j-1
            if (data(indx(j)).le.datap) then
               exit Q4
            endif
         enddo Q4
         if (i.lt.j) then
            indext=indx(i)
            indx(i)=indx(j)
            indx(j)=indext
            cycle Q3
         else
            if (r-j .ge. i-l .and. i-l .gt. m) then
               istk=istk+1
               lstk(istk)=j+1
               rstk(istk)=r
               r=i-1
            else if (i-l .gt. r-j .and. r-j .gt. m) then
               istk=istk+1
               lstk(istk)=l
               rstk(istk)=i-1
               l=j+1
            else if (r-j .gt. m) then
               l=j+1
            else if (i-l .gt. m) then
               r=i-1
            else
               if (istk.lt.1) then
                  exit TOP
               endif
               l=lstk(istk)
               r=rstk(istk)
               istk=istk-1
            endif
            cycle TOP
         endif
      enddo Q3
      exit TOP
   enddo TOP
endif
do i=2,n
   if (data(indx(i-1)) .gt. data(indx(i))) then
      indexp=indx(i)
      datap=data(indexp)
      p=i-1
      INNER: do
         indx(p+1) = indx(p)
         p=p-1
         if (p.le.0)then
            exit INNER
         endif
         if (data(indx(p)).le.datap)then
            exit INNER
         endif
      enddo INNER
      indx(p+1) = indexp
   endif
enddo
end subroutine sort_quick_rx_${TYPE}_${KIND}_${IKIND}
$ENDPARCEL 
$!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$parcel sort_ikind
$undefine character

$SET TYPE integer 
$SET KIND int8 
$POST SORT_QUICK_RX 
$SET KIND int16 
$POST SORT_QUICK_RX 
$SET KIND int32 
$POST SORT_QUICK_RX 
$SET KIND int64 
$POST SORT_QUICK_RX 
$SET TYPE real 
$SET KIND real32 
$POST SORT_QUICK_RX 
$SET KIND real64 
$POST SORT_QUICK_RX 
$define character 
$SET TYPE character 
$SET KIND ascii 
$POST SORT_QUICK_RX 

$undefine character
$endparcel
$!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
$set ikind int32
$post sort_ikind
$set ikind int64
$post sort_ikind

end module M_sort

program demo_sort_quick_rx
use M_sort, only : sort_quick_rx
implicit none
integer,parameter            :: isz=10000
real                         :: rr(isz)
integer                      :: ii(isz)
integer                      :: i
   write(*,*)'initializing array with ',isz,' random numbers'
   CALL RANDOM_NUMBER(RR)
   rr=rr*450000.0
   write(*,*)'sort real array with sort_quick_rx(3f)'
   call sort_quick_rx(rr,ii)
   write(*,*)'checking index of sort_quick_rx(3f)'
   do i=1,isz-1
      if(rr(ii(i)).gt.rr(ii(i+1)))then
         write(*,*)'Error in sorting reals small to large ', &
         & i,rr(ii(i)),rr(ii(i+1))
      endif
   enddo
   write(*,*)'test of sort_quick_rx(3f) complete'
   ! use the index array to actually move the input array into a sorted
   ! order
   rr=rr(ii)
   do i=1,isz-1
      if(rr(i).gt.rr(i+1))then
         write(*,*)'Error in sorting reals small to large ', &
         & i,rr(i),rr(i+1)
      endif
   enddo
   write(*,*)'test of sort_quick_rx(3f) complete'

end program demo_sort_quick_rx