module sort

  use :: iso_c_binding

  implicit none

  private
  public :: Qsort

  integer, parameter :: dp = c_double
  integer, parameter :: i4 = c_int
  integer, parameter :: i8 = c_int64_t

  type Limits
      integer :: Ileft, Iright
  end type Limits

  interface Qsort
    module procedure :: Qsort_i4
    module procedure :: Qsort_i8
    module procedure :: Qsort_double
  end interface Qsort

  interface ChoosePiv
    module procedure :: ChoosePiv_i4
    module procedure :: ChoosePiv_i8
    module procedure :: ChoosePiv_double
  end interface ChoosePiv

  interface InsrtLC
    module procedure :: InsrtLC_i4
    module procedure :: InsrtLC_i8
    module procedure :: InsrtLC_double
  end interface InsrtLC

  interface Partition
    module procedure :: Partition_i4
    module procedure :: Partition_i8
    module procedure :: Partition_double
  end interface Partition

  interface Swap
    module procedure :: Swap_i4
    module procedure :: Swap_i8
    module procedure :: Swap_double
  end interface Swap

  interface Reverse
    module procedure :: Reverse_i4
    module procedure :: Reverse_i8
    module procedure :: Reverse_double
  end interface Reverse

contains

  ! ***********************************
  ! * Sort Array X(:) in ascendent order 
  ! * If present Ipt, a pointer with the 
  ! * changes is returned in Ipt.
  ! * If present order, +1 gives ascendent order (default)
  ! *                   -1 gives descendent order
  ! ***********************************
  subroutine Qsort_i4(X, Ipt, order)
    implicit none
    integer(kind=c_int), dimension(:), intent(inout) :: X
    integer(kind=c_int), dimension(:), intent(out), optional :: Ipt(size(X))
    integer(kind=c_int), intent(in), optional :: order
  
    ! For a list with Isw number of elements or
    ! less use Insrt
    integer, parameter :: Isw = 10
      
    integer :: I, Ipvn, Ileft, Iright, ISpos, ISmax
    integer, dimension(:), allocatable :: IIpt
    type(Limits), dimension(:), allocatable :: Stack

    integer :: N

    N = size(X)
      
    allocate(Stack(N))

    Stack(:)%Ileft = 0
    if (present(Ipt)) then
      Ipt = (/ (I, I=1,N) /)

      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
        
      do while (Stack(ISpos)%Ileft /= 0)

          Ileft = Stack(ISPos)%Ileft
          Iright = Stack(ISPos)%Iright
          if (Iright-Ileft <= Isw) Then
              call InsrtLC(X, Ipt, Ileft,Iright)
              ISpos = ISPos + 1
          else
              Ipvn = ChoosePiv(X, Ileft, Iright)
              Ipvn = Partition(X, Ileft, Iright, Ipvn, Ipt)
            
              Stack(ISmax+1)%Ileft = Ileft
              Stack(ISmax+1) %Iright = Ipvn-1
              Stack(ISmax+2)%Ileft = Ipvn + 1
              Stack(ISmax+2)%Iright = Iright
              ISpos = ISpos + 1
              ISmax = ISmax + 2
          end if
      end do

    else
      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
      
      allocate(IIpt(10))
      do while (Stack(ISpos)%Ileft /= 0)
          Ileft = Stack(ISPos)%Ileft
          Iright = Stack(ISPos)%Iright
          if (Iright-Ileft <= Isw) Then
              call InsrtLC(X, IIpt, Ileft, Iright)
              ISpos = ISPos + 1
          else
              Ipvn = ChoosePiv(X, Ileft, Iright)
              Ipvn = Partition(X, Ileft, Iright, Ipvn)
              
              Stack(ISmax+1)%Ileft = Ileft
              Stack(ISmax+1) %Iright = Ipvn-1
              Stack(ISmax+2)%Ileft = Ipvn + 1
              Stack(ISmax+2)%Iright = Iright
              ISpos = ISpos + 1
              ISmax = ISmax + 2
          end if
      end do
      deallocate(IIpt)

    end if

    deallocate(Stack)

    if (present(order)) then
      if (order==-1) then
        call Reverse(X)
        call Reverse(Ipt)
      end if
    end if
  
    return

  end subroutine Qsort_i4

  subroutine Qsort_i8(X, Ipt, order)
    implicit none
    integer(kind=c_int64_t), dimension(:), intent(inout) :: X
    integer(kind=c_int64_t), dimension(:), intent(out), optional :: Ipt(size(X))
    integer(kind=c_int), intent(in), optional :: order
  
    ! For a list with Isw number of elements or
    ! less use Insrt
    integer, parameter :: Isw = 10
      
    integer :: I, Ipvn, Ileft, Iright, ISpos, ISmax
    integer(kind=c_int64_t), dimension(:), allocatable :: IIpt
    type(Limits), dimension(:), allocatable :: Stack

    integer :: N

    N = size(X)
      
    allocate(Stack(N))

    Stack(:)%Ileft = 0
    if (present(Ipt)) then
      Ipt = (/ (I, I=1,N) /)

      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
        
      do while (Stack(ISpos)%Ileft /= 0)

        Ileft = Stack(ISPos)%Ileft
        Iright = Stack(ISPos)%Iright
        if (Iright-Ileft <= Isw) Then
          call InsrtLC(X, Ipt, Ileft,Iright)
          ISpos = ISPos + 1
        else
          Ipvn = ChoosePiv(X, Ileft, Iright)
          Ipvn = Partition(X, Ileft, Iright, Ipvn, Ipt)
        
          Stack(ISmax+1)%Ileft = Ileft
          Stack(ISmax+1) %Iright = Ipvn-1
          Stack(ISmax+2)%Ileft = Ipvn + 1
          Stack(ISmax+2)%Iright = Iright
          ISpos = ISpos + 1
          ISmax = ISmax + 2
        end if
      end do

    else
      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
      
      allocate(IIpt(10))
      do while (Stack(ISpos)%Ileft /= 0)
        Ileft = Stack(ISPos)%Ileft
        Iright = Stack(ISPos)%Iright
        if (Iright-Ileft <= Isw) Then
          call InsrtLC(X, IIpt, Ileft, Iright)
          ISpos = ISPos + 1
        else
          Ipvn = ChoosePiv(X, Ileft, Iright)
          Ipvn = Partition(X, Ileft, Iright, Ipvn)
            
          Stack(ISmax+1)%Ileft = Ileft
          Stack(ISmax+1) %Iright = Ipvn-1
          Stack(ISmax+2)%Ileft = Ipvn + 1
          Stack(ISmax+2)%Iright = Iright
          ISpos = ISpos + 1
          ISmax = ISmax + 2
        end if
      end do
      deallocate(IIpt)

    end if

    deallocate(Stack)

    if (present(order)) then
      if (order==-1) then
        call Reverse(X)
        call Reverse(Ipt)
      end if
    end if
  
    return

  end subroutine Qsort_i8

  subroutine Qsort_double(X, Ipt, order)
    implicit none
    real(kind=c_double), dimension(:), intent(inout) :: X
    integer(kind=c_int), dimension(:), intent(out), optional :: Ipt(size(X))
    integer(kind=c_int), intent(in), optional :: order
  
    ! For a list with Isw number of elements or
    ! less use Insrt
    integer, parameter :: Isw = 10
      
    integer :: I, Ipvn, Ileft, Iright, ISpos, ISmax
    integer, dimension(:), allocatable :: IIpt
    type(Limits), dimension(:), allocatable :: Stack

    integer :: N

    N = size(X)
      
    allocate(Stack(N))

    Stack(:)%Ileft = 0
    if (present(Ipt)) then
      Ipt = (/ (I, I=1,N) /)

      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
        
      do while (Stack(ISpos)%Ileft /= 0)

        Ileft = Stack(ISPos)%Ileft
        Iright = Stack(ISPos)%Iright
        if (Iright-Ileft <= Isw) Then
          call InsrtLC(X, Ipt, Ileft,Iright)
          ISpos = ISPos + 1
        else
          Ipvn = ChoosePiv(X, Ileft, Iright)
          Ipvn = Partition(X, Ileft, Iright, Ipvn, Ipt)
        
          Stack(ISmax+1)%Ileft = Ileft
          Stack(ISmax+1) %Iright = Ipvn-1
          Stack(ISmax+2)%Ileft = Ipvn + 1
          Stack(ISmax+2)%Iright = Iright
          ISpos = ISpos + 1
          ISmax = ISmax + 2
        end if
      end do

    else
      ! Iniitialize the stack
      Ispos = 1
      Ismax = 1
      Stack(ISpos)%Ileft  = 1
      Stack(ISpos)%Iright = N
      
      allocate(IIpt(10))
      do while (Stack(ISpos)%Ileft /= 0)
        Ileft = Stack(ISPos)%Ileft
        Iright = Stack(ISPos)%Iright
        if (Iright-Ileft <= Isw) Then
          call InsrtLC(X, IIpt, Ileft, Iright)
          ISpos = ISPos + 1
        else
          Ipvn = ChoosePiv(X, Ileft, Iright)
          Ipvn = Partition(X, Ileft, Iright, Ipvn)
            
          Stack(ISmax+1)%Ileft = Ileft
          Stack(ISmax+1) %Iright = Ipvn-1
          Stack(ISmax+2)%Ileft = Ipvn + 1
          Stack(ISmax+2)%Iright = Iright
          ISpos = ISpos + 1
          ISmax = ISmax + 2
        end if
      end do
      deallocate(IIpt)

    end if

    deallocate(Stack)

    if (present(order)) then
      if (order==-1) then
        call Reverse(X)
        call Reverse(Ipt)
      end if
    end if
  
    return

  end subroutine Qsort_double
        
  ! ***********************************
  ! * Choose a Pivot element from XX(Ileft:Iright)
  ! * for Qsort. This routine chooses the median
  ! * of the first, last and mid element of the 
  ! * list.
  ! ***********************************
  function ChoosePiv_i4(XX, IIleft, IIright) result(IIpv)
    
    implicit none
    integer :: IIpv        
    integer(kind=c_int), dimension(:), intent(in) :: XX
    integer, intent(in) :: IIleft, IIright
          
    integer(kind=c_int) :: XXcp(3)
    integer :: IIpt(3), IImd
          
    IImd = Int((IIleft+IIright)/2)
    XXcp(1) = XX(IIleft)
    XXcp(2) = XX(IImd)
    XXcp(3) = XX(IIright)
    IIpt = (/1,2,3/)
    
    call InsrtLC(XXcp, IIpt, 1, 3)
          
    select case (IIpt(2))
      case (1)
        IIpv = IIleft
      case (2)
        IIpv = IImd
      case (3)
        IIpv = IIright
    end select
    
    return
  end function ChoosePiv_i4

  function ChoosePiv_i8(XX, IIleft, IIright) result(IIpv)
    
    implicit none
    integer :: IIpv        
    integer(kind=c_int64_t), dimension(:), intent(in) :: XX
    integer, intent(in) :: IIleft, IIright
          
    integer(kind=c_int64_t) :: XXcp(3)
    integer :: IImd
    integer(kind=c_int64_t) :: IIpt(3)
          
    IImd = Int((IIleft+IIright)/2)
    XXcp(1) = XX(IIleft)
    XXcp(2) = XX(IImd)
    XXcp(3) = XX(IIright)
    IIpt = (/1,2,3/)
    
    call InsrtLC(XXcp, IIpt, 1, 3)
          
    select case (IIpt(2))
      case (1)
        IIpv = IIleft
      case (2)
        IIpv = IImd
      case (3)
        IIpv = IIright
    end select
    
    return
  end function ChoosePiv_i8

  function ChoosePiv_double(XX, IIleft, IIright) result(IIpv)
    
    implicit none
    integer :: IIpv        
    real(kind=c_double), dimension(:), intent(in) :: XX
    integer, intent(in) :: IIleft, IIright
          
    real(kind=c_double) :: XXcp(3)
    integer :: IIpt(3), IImd
          
    IImd = Int((IIleft+IIright)/2)
    XXcp(1) = XX(IIleft)
    XXcp(2) = XX(IImd)
    XXcp(3) = XX(IIright)
    IIpt = (/1,2,3/)
    
    call InsrtLC(XXcp, IIpt, 1, 3)
          
    select case (IIpt(2))
      case (1)
        IIpv = IIleft
      case (2)
        IIpv = IImd
      case (3)
        IIpv = IIright
    end select
    
    return
  end function ChoosePiv_double
    
  ! ***********************************
  ! * Perform an insertion sort of the list 
  ! * XX(:) between index values IIl and IIr.
  ! * IIpt(:) returns the permutations
  ! * made to sort.
  ! ***********************************
  subroutine InsrtLC_i4(XX, IIpt, IIl, IIr)
    implicit none    
    integer(kind=c_int), dimension(:), intent(inout) :: XX
    integer, dimension(:), intent(inout) :: IIpt
    integer, intent(in) :: IIl, IIr
          
    integer(kind=c_int) :: RRtmp
    integer :: II, JJ
    
    do II = IIl+1, IIr
      RRtmp = XX(II)
      do JJ = II-1, 1, -1
        if (RRtmp < XX(JJ)) then
          XX(JJ+1) = XX(JJ)
          call Swap(IIpt, JJ, JJ+1)
        else
          exit
        end if
      end do
      XX(JJ+1) = RRtmp
    end do
          
    return
  end subroutine InsrtLC_i4

  subroutine InsrtLC_i8(XX, IIpt, IIl, IIr)
    implicit none    
    integer(kind=c_int64_t), dimension(:), intent(inout) :: XX
    integer(kind=c_int64_t), dimension(:), intent(inout) :: IIpt
    integer, intent(in) :: IIl, IIr
          
    integer(kind=c_int64_t) :: RRtmp
    integer :: II, JJ
    
    do II = IIl+1, IIr
      RRtmp = XX(II)
      do JJ = II-1, 1, -1
        if (RRtmp < XX(JJ)) then
          XX(JJ+1) = XX(JJ)
          call Swap(IIpt, JJ, JJ+1)
        else
          exit
        end if
      end do
      XX(JJ+1) = RRtmp
    end do
          
    return
  end subroutine InsrtLC_i8

  subroutine InsrtLC_double(XX, IIpt, IIl, IIr)
    implicit none    
    real(kind=c_double), dimension(:), intent(inout) :: XX
    integer, dimension(:), intent(inout) :: IIpt
    integer, intent(in) :: IIl, IIr
          
    real(kind=c_double) :: RRtmp
    integer :: II, JJ
    
    do II = IIl+1, IIr
      RRtmp = XX(II)
      do JJ = II-1, 1, -1
        if (RRtmp < XX(JJ)) then
          XX(JJ+1) = XX(JJ)
          call Swap(IIpt, JJ, JJ+1)
        else
          exit
        end if
      end do
      XX(JJ+1) = RRtmp
    end do
          
    return
  end subroutine InsrtLC_double
    
  ! ***********************************
  ! * This routine arranges the array X
  ! * between the index values Ileft and Iright
  ! * positioning elements smallers than
  ! * X(Ipv) at the left and the others 
  ! * at the right.
  ! * Internal routine used by Qsort.
  ! ***********************************
  function Partition_i4(X, Ileft, Iright, Ipv, Ipt) Result (Ipvfn)
    implicit none    
    integer(kind=c_int), dimension(:), intent(inout) :: X
    integer, intent(in) :: Ileft, Iright, Ipv
    integer, intent(inout), dimension(:), optional :: Ipt
    integer :: Ipvfn

    integer(kind=c_int) :: Rpv
    integer :: I

    Rpv = X(Ipv)
    call Swap(X, Ipv, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipv, Iright)
    Ipvfn = Ileft

    if (present(Ipt))  then
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          call Swap(Ipt, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    else
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    end if

    call Swap(X, Ipvfn, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipvfn, Iright)

    return
  end function Partition_i4

  function Partition_i8(X, Ileft, Iright, Ipv, Ipt) Result (Ipvfn)
    implicit none    
    integer(kind=c_int64_t), dimension(:), intent(inout) :: X
    integer, intent(in) :: Ileft, Iright, Ipv
    integer(kind=c_int64_t), intent(inout), dimension(:), optional :: Ipt
    integer :: Ipvfn

    integer(kind=c_int64_t) :: Rpv
    integer :: I

    Rpv = X(Ipv)
    call Swap(X, Ipv, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipv, Iright)
    Ipvfn = Ileft

    if (present(Ipt))  then
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          call Swap(Ipt, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    else
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    end if

    call Swap(X, Ipvfn, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipvfn, Iright)

    return
  end function Partition_i8

  function Partition_double(X, Ileft, Iright, Ipv, Ipt) Result (Ipvfn)
    implicit none    
    real(kind=c_double), dimension(:), intent(inout) :: X
    integer, intent(in) :: Ileft, Iright, Ipv
    integer, intent(inout), dimension(:), optional :: Ipt
    integer :: Ipvfn

    real(kind=c_double) :: Rpv
    integer :: I

    Rpv = X(Ipv)
    call Swap(X, Ipv, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipv, Iright)
    Ipvfn = Ileft

    if (present(Ipt))  then
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          call Swap(Ipt, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    else
      do I = Ileft, Iright-1
        if (X(I) <= Rpv) then
          call Swap(X, I, Ipvfn)
          Ipvfn = Ipvfn + 1
        end if
      end do
    end if

    call Swap(X, Ipvfn, Iright)
    if (present(Ipt)) call Swap(Ipt, Ipvfn, Iright)

    return
  end function Partition_double
    
  ! ***********************************
  ! * Swaps elements I and J of array X(:). 
  ! ***********************************
  subroutine Swap_i4(X, I, J)
    implicit none
    integer(kind=c_int), dimension(:), intent(inout) :: X
    integer, intent(in) :: I, J

    integer(kind=c_int) :: Itmp

    Itmp = X(I)
    X(I) = X(J)
    X(J) = Itmp

    return
  end subroutine Swap_i4

  subroutine Swap_i8(X, I, J)
    implicit none
    integer(kind=c_int64_t), dimension(:), intent(inout) :: X
    integer, intent(in) :: I, J

    integer(kind=c_int64_t) :: Itmp

    Itmp = X(I)
    X(I) = X(J)
    X(J) = Itmp

    return
  end subroutine Swap_i8

  subroutine Swap_double(X, I, J)
    implicit none
    real(kind=c_double), dimension(:), intent(inout) :: X
    integer, intent(in) :: I, J

    real(kind=c_double) :: Itmp

    Itmp = X(I)
    X(I) = X(J)
    X(J) = Itmp

    return
  end subroutine Swap_double

  subroutine Reverse_i4(X)
    implicit none
    integer(kind=c_int), dimension(:), intent(inout) :: X

    X = X(ubound(X, dim=1):lbound(X, dim=1):-1)
    return
  end subroutine Reverse_i4

  subroutine Reverse_i8(X)
    implicit none
    integer(kind=c_int64_t), dimension(:), intent(inout) :: X

    X = X(ubound(X, dim=1):lbound(X, dim=1):-1)
    return
  end subroutine Reverse_i8

  subroutine Reverse_double(X)
    implicit none
    real(kind=c_double), dimension(:), intent(inout) :: X

    X = X(ubound(X, dim=1):lbound(X, dim=1):-1)
    return
  end subroutine Reverse_double

end module sort