! sortmod.f90
! module procedure written by Tim Mitchell on 11.7.00 and thereafter
! QuickSort etc
! BubbleSort, MedianValue, RankVector are in old versions

module SortMod

implicit none

contains

!*******************************************************************************
! quick sort text in ascending order - call this procedure
! call with Text, returns in Position the alphabetical order of the indices of Text

subroutine QuickSortText (Text,Position,Restricted)

real, pointer, dimension (:) 			:: SubSet
real, allocatable, dimension (:,:) 		:: AsciiCode

integer, pointer, dimension (:) 		:: Position
integer, allocatable, dimension (:) 		:: Order

character (len=80), pointer, dimension (:) 	:: Text

integer, intent(in), optional 			:: Restricted	! contains the max char of Text to consider

real, parameter :: MissVal = -999.0

integer :: AllocStat
integer :: XDatum, XChar, XSubSet, XPlace
integer :: DataN,  CharN, SubSetN
integer :: MaxPlace, CurrentAsciiCode,CurrentDatum,CurrentPlace, LastAsciiCode,PlacesToAdd

character (len=80) :: DatumText

DataN = size (Text,1)

if (present(Restricted)) then
  CharN = Restricted
else
  CharN = 80
end if

allocate (AsciiCode  (DataN,CharN), &
	  Order      (DataN),       &
	  Position   (DataN),       stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSortText: Allocation failure: most #####"
Order = 1

do XDatum = 1, DataN  
  DatumText = trim(adjustl(Text(XDatum)))
  do XChar = 1, CharN
    AsciiCode (XDatum,XChar) = real(iachar(DatumText(XChar:XChar))) + (real(XDatum)/10000.0)
  end do
end do

do XChar = 1, CharN
 do XPlace = DataN, 1, -1
  SubSetN = count (Order.EQ.XPlace)
  
  if (SubSetN.GT.1) then  
   allocate (SubSet (SubSetN),stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSortText: Allocation failure: SubSet #####"

   XSubSet = 0
   do XDatum = 1, DataN
     if (Order(XDatum).EQ.XPlace) then
    	XSubSet=XSubSet+1 ; SubSet(XSubSet) = AsciiCode(XDatum,XChar)
     end if
   end do
   
   call REFSOR (SubSet)
   
   CurrentPlace = XPlace-1 ; LastAsciiCode = -1 ; PlacesToAdd = 1
   do XSubSet = 1, SubSetN
     CurrentAsciiCode = floor (SubSet(XSubSet))
     CurrentDatum     = nint ((SubSet(XSubSet)-real(CurrentAsciiCode))*10000.0)
     
     if (CurrentAsciiCode.GT.LastAsciiCode) then
     	CurrentPlace  = CurrentPlace + PlacesToAdd
     	PlacesToAdd   = 1
     	LastAsciiCode = CurrentAsciiCode
     else
        PlacesToAdd  = PlacesToAdd + 1
     end if
     
     Order(CurrentDatum) = CurrentPlace
   end do
     
   deallocate (SubSet,stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSortText: Deallocation failure: Column #####"
  end if
 end do
end do

do XDatum = 1, DataN  
  Position(Order(XDatum)) = XDatum
end do

deallocate (AsciiCode,Order,stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSortText: Deallocation failure: AsciiCode #####"

end subroutine QuickSortText 

!*******************************************************************************
! quick sort in ascending order - call this procedure
! if Order or OrderValid is set, Reals/Ints are left unchanged
! the difference between Order and OrderValid is that
!	Order=intarr(DataN) and is declared here and passed to calling program
!	OrderValid=intarr(at least ValidN) and is already declared

subroutine QuickSort (Reals,Ints,Order,OrderValid,NMiss)

logical, pointer, dimension (:)			:: Found
real, pointer, dimension (:)			:: Data,ValidData
real, pointer, dimension (:), optional 		:: Reals
integer, pointer, dimension (:), optional 	:: Order,OrderValid,Ints

integer, intent (out), optional		:: NMiss

real, parameter :: MissVal = -999.0
integer :: XDatum, XValid, ValidN, XSort, AllocStat, DataN, MissN

if (present(Reals)) then
  DataN = size(Reals)
  allocate (Data (DataN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Allocation failure: Data #####"
  do XDatum = 1, DataN
    Data(XDatum) = Reals(XDatum)
  end do
else if (present(Ints)) then
  DataN = size(Ints)
  allocate (Data (DataN), stat=AllocStat)
  if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Allocation failure: Data #####"
  do XDatum = 1, DataN
    Data(XDatum) = real(Ints(XDatum))
  end do
else
  print*, "  > ##### ERROR: QuickSort: no data supplied! #####"
end if

MissN = 0
do XDatum = 1, DataN
  if (Data(XDatum).EQ.MissVal) MissN = MissN + 1
end do
ValidN = DataN - MissN
if (present(NMiss)) NMiss = MissN

allocate (ValidData (ValidN), stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Allocation failure #####"

XValid = 0
do XDatum = 1, DataN
  if (Data(XDatum).NE.MissVal) then
  	XValid = XValid + 1
	ValidData(XValid) = Data(XDatum)
  end if
end do

if (ValidN.GT.1) call REFSOR (ValidData)

if (present(Order).OR.present(OrderValid)) then
 allocate (Found(DataN), stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Allocation failure #####"
 Found=.FALSE.
 
 if (present(Order)) then
   allocate (Order(DataN), stat=AllocStat)
   if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Allocation failure #####"
   Order = MissVal

   do XSort = 1, ValidN
    XDatum = 0
    do
      XDatum = XDatum + 1
      if (ValidData(XSort).EQ.Data(XDatum).AND.(.NOT.Found(XDatum))) then
        Found(XDatum) = .TRUE.
        Order(XSort)  = XDatum
      end if
      if (Order(XSort).NE.MissVal) exit
    end do
   end do
 else
   OrderValid = nint(MissVal)
   if (size(OrderValid).LT.ValidN) print*, "  > ##### ERROR: QuickSort: OrderValid #####"
   
   do XSort = 1, ValidN
    XDatum = 0
    do
      XDatum = XDatum + 1
      if (ValidData(XSort).EQ.Data(XDatum).AND.(.NOT.Found(XDatum))) then
        Found(XDatum) = .TRUE.
        OrderValid(XSort) = XDatum
      end if
      if (OrderValid(XSort).NE.nint(MissVal)) exit
    end do
   end do
 end if
 
 deallocate (Found, stat=AllocStat)
 if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Dealloc: Found #####"
else
 if (present(Reals)) then 
   do XDatum = 1, ValidN
     Reals (XDatum) = ValidData (XDatum)
   end do
   do XDatum = (ValidN+1), DataN
     Reals (XDatum) = MissVal
   end do
 else if (present(Ints)) then
   do XDatum = 1, ValidN
     Ints (XDatum) = ValidData (XDatum)
   end do
   do XDatum = (ValidN+1), DataN
     Ints (XDatum) = MissVal
   end do
 end if
end if

deallocate (ValidData,Data, stat=AllocStat)
if (AllocStat.NE.0) print*, "  > ##### ERROR: QuickSort: Deallocation failure #####"

end subroutine QuickSort

!*******************************************************************************
! quick sort in ascending order - the actual procedure - do not call directly

Subroutine REFSOR (XDONT)
!  Sorts XDONT into ascending order - Quicksort
! __________________________________________________________
!  Quicksort chooses a "pivot" in the set, and explores the
!  array from both ends, looking for a value > pivot with the
!  increasing index, for a value <= pivot with the decreasing
!  index, and swapping them when it has found one of each. 
!  The array is then subdivided in 2 ([3]) subsets:
!  { values <= pivot} {pivot} {values > pivot}
!  One then call recursively the program to sort each subset.
!  When the size of the subarray is small enough, one uses an
!  insertion sort that is faster for very small sets.
!  Michel Olagnon - Apr. 2000
! __________________________________________________________
!
!      Real, Dimension (:), Intent (InOut) :: XDONT
!  modified by Tim Mitchell, 11.7.00
      real, dimension (:), pointer :: XDONT
! __________________________________________________________
      Integer :: IDEB, IFIN
!
      IDEB = 1
      IFIN = Size (XDONT)
      Call l_subsor (IDEB, IFIN)
      Call l_inssor
      Return
Contains      
   Recursive Subroutine l_subsor (IDEB1, IFIN1)
!  Sorts XDONT from IDEB1 to IFIN1
! __________________________________________________________
      Integer, Intent (In) :: IDEB1, IFIN1
! __________________________________________________________
      Integer, Parameter :: NINS = 16 ! Max for insertion sort 
      Integer :: ICRS, IDEB, IDCR, IFIN, IMIL
      Real (Kind(XDONT)) :: XPIV, XWRK
!
      IDEB = IDEB1
      IFIN = IFIN1
!
!  If we don't have enough values to make it worth while, we leave
!  them unsorted, and the final insertion sort will take care of them
!
      If ((IFIN - IDEB) > NINS) Then
         IMIL = (IDEB+IFIN) / 2
!
!  One chooses a pivot, median of 1st, last, and middle values
!
         If (XDONT(IMIL) < XDONT(IDEB)) Then
            XWRK = XDONT (IDEB)
            XDONT (IDEB) = XDONT (IMIL)
            XDONT (IMIL) = XWRK
         End If
         If (XDONT(IMIL) > XDONT(IFIN)) Then
            XWRK = XDONT (IFIN)
            XDONT (IFIN) = XDONT (IMIL)
            XDONT (IMIL) = XWRK
            If (XDONT(IMIL) < XDONT(IDEB)) Then
               XWRK = XDONT (IDEB)
               XDONT (IDEB) = XDONT (IMIL)
               XDONT (IMIL) = XWRK
            End If
         End If
         XPIV = XDONT (IMIL)
!
!  One exchanges values to put those > pivot in the end and
!  those <= pivot at the beginning
!
         ICRS = IDEB
         IDCR = IFIN
         ECH2: Do
            Do
               ICRS = ICRS + 1
               If (ICRS >= IDCR) Then
!
!  the first  >  pivot is IDCR
!  the last   <= pivot is ICRS-1
!  Note: If one arrives here on the first iteration, then
!        the pivot is the maximum of the set, the last value is equal
!        to it, and one can reduce by one the size of the set to process,
!        as if XDONT (IFIN) > XPIV
!
                  Exit ECH2
!
               End If
               If (XDONT(ICRS) > XPIV) Exit
            End Do
            Do
               If (XDONT(IDCR) <= XPIV) Exit
               IDCR = IDCR - 1
               If (ICRS >= IDCR) Then
!
!  The last value < pivot is always ICRS-1
!
                  Exit ECH2
               End If
            End Do
!
            XWRK = XDONT (IDCR)
            XDONT (IDCR) = XDONT (ICRS)
            XDONT (ICRS) = XWRK
         End Do ECH2
!
!  One now sorts each of the two sub-intervals
!
         Call l_subsor (IDEB1, ICRS-1)
         Call l_subsor (IDCR, IFIN1)
      End If
      Return
   End Subroutine l_subsor
   Subroutine l_inssor
!  Sorts XDONT into increasing order (Insertion sort)
!  This subroutine uses insertion sort. It does not use any
!  work array and is faster when XDONT is of very small size
!  (< 20), or already almost sorted, so it is used in a final
!  pass when the partial quicksorting has left a sequence
!  of small subsets and that sorting is only necessary within
!  each subset to complete the process.
!  Michel Olagnon - Apr. 2000
! __________________________________________________________
      Integer :: ICRS, IDCR
      Real (Kind(XDONT)) :: XWRK
!
      Do ICRS = 2, Size (XDONT)
         XWRK = XDONT (ICRS)
         If (XWRK >= XDONT(ICRS-1)) Cycle
         XDONT (ICRS) = XDONT (ICRS-1)
         Do IDCR = ICRS - 2, 1, - 1
            If (XWRK >= XDONT(IDCR)) Exit
            XDONT (IDCR+1) = XDONT (IDCR)
         End Do
         XDONT (IDCR+1) = XWRK
      End Do
!
      Return
!
   End Subroutine l_inssor
End Subroutine REFSOR

!*******************************************************************************

end module SortMod