! pair_sort.f90
!
! always  Order n lg n  sort, reliably fast
! reference: P 147 of Introduction to Algorithms by Cormen et. al.
!
! 11/5/95 JSS converted from Ada with minimum editing
! define dervied type and sort on random numbers

module pair_sort
  implicit none
  type deck
    integer :: card
    real :: random
  end type deck

  type (deck), dimension(0:19) :: cards

contains
  subroutine PairInit
    integer :: i
    integer :: j
    integer, dimension(1) :: seed

    j = -1
    do i=0,19
      if (modulo(i,2)==0) j=j+1
      cards(i)%card = j
    end do
    call system_clock(seed(1))
    call random_seed(put=seed)
  end subroutine PairInit

  subroutine PairDeal
    integer :: i

    do i=0,19
      call random_number(cards(i)%random) 
    end do
    call PairSort(cards)
  end subroutine PairDeal

  subroutine PairSort ( A )
    implicit none
    type (deck), dimension(:), intent(inout) :: A
    type (deck) :: T
    integer :: Heap_Size
    integer :: I

    if (size(A) <= 1) then
      return
    end if
    call Build_Heap(A)
    do I=ubound(A,1), lbound(A,1), -1
      T = A(lbound(A,1));
      A(lbound(A,1)) = A(I)
      A(I) = T
      Heap_Size = Heap_Size-1
      call Heapify(A,lbound(A,1))
    end do

  contains
    recursive subroutine Heapify(A, I)
      type (deck), dimension(:), intent(inout) :: A
      integer, intent(in) :: I
      type (deck) :: T
      integer :: Largest
      integer :: L
      integer :: R

      L = lbound(A,1) + 2*(I-lbound(A,1))
      R = L+1
      if (L <= Heap_Size .and. A(L)%random > A(I)%random) then
        Largest = L
      else
        Largest = I
      end if
      if (R <= Heap_Size .and. A(R)%random > A(Largest)%random) then
        Largest = R
      end if
      if (Largest /= I) then
        T = A(I)
        A(I) = A(Largest)
        A(Largest) = T
        call Heapify(A, Largest)
      end if
    end subroutine Heapify

    subroutine Build_Heap ( A )
      type (deck), dimension(:), intent(inout) :: A
      integer :: I

      Heap_Size = size(A,1)
      do I=lbound(A,1)+size(A,1)/2 ,lbound(A,1) ,-1
        call Heapify(A, I);
      end do
    end subroutine Build_Heap;
  end subroutine PairSort;
end module Pair_Sort