!     Last change:  DOS   3 Aug 2000    6:08 pm
!                                       *** copyright 2000 ***
! *** filename hsort.f95                *** John F. Monahan **
!                                       **********************
      program phsort
!  TEST OF HEAPSORT ALGORITHM HSORT
      implicit none
      real, dimension(14400) :: x
      real ran
      integer i,j,k,n,iflag
!
  20  FORMAT(//' Test Problems with Sample Size',i8)
  21  FORMAT(/' Simple Case -- All Positive, N=',i8,'  OK if zero',i2)
  23  FORMAT(/' Simple Case -- All Negative, N=',i8,'  OK if zero',i2)
  25  FORMAT(/' Simple Case -- Mixed Signs, N=',i8,'  OK if zero',i2)
  27  FORMAT(/' Mixed Signs and Some Ties, N=',i8,'  OK if zero',i2)
  22  FORMAT(2X,8F9.4)
  26  FORMAT(2X,8F9.6)
  28  FORMAT(2X,12F6.2)
!                   open output file
      open( unit=6, file='hsort.out' )
!                   initialize uniform generator
      x(1) = ran(5151917)
!                   let's check out these hundred and twenty cases
      do i = 1,120
      n = i 
      if( i .gt. 10 ) n = i*i
      write(6,20) n
!                   first all positive
      do j = 1,n
      x(j) = 5.*ran(j)
      end do  ! loop on j
!                   call sorting algorithm and write out results
      call hsort(x,n)
!                   test
      iflag = 0
      do j = 2,n
      if( x(j-1) .gt. x(j) ) iflag = 1
      end do  ! loop on j
      write(6,21) n,iflag
      if( n .lt. 20 ) write(6,22) (x(j),j=1,n)
!                   now all negative
      do j = 1,n
      x(j) = -100.*ran(j)
      end do  ! loop on j
!                   call sorting algorithm and write out results
      call hsort(x,n)
!                   test
      iflag = 0
      do j = 2,n
      if( x(j-1) .gt. x(j) ) iflag = 1
      end do  ! loop on j
      write(6,23) n,iflag
      if( n .lt. 20 ) write(6,22) (x(j),j=1,n)
!                   now mixed positive and negative and smaller
      do j = 1,n
      x(j) = (ran(j)-.3)/100.
      end do  ! loop on j
!                   call sorting algorithm and write out results
      call hsort(x,n)
!                   test
      iflag = 0
      do j = 2,n
      if( x(j-1) .gt. x(j) ) iflag = 1
      end do  ! loop on j
      write(6,25) n,iflag
      if( n .lt. 20 ) write(6,26) (x(j),j=1,n)
!                   now introduce some ties
      do j = 1,n
      k = int( 8.*(ran(j)-.3) )
      x(j) = real(k)/10.
      end do  ! loop on j
!                   call sorting algorithm and write out results
      call hsort(x,n)
!                   test
      iflag = 0
      do j = 2,n
      if( x(j-1) .gt. x(j) ) iflag = 1
      end do  ! loop on j
      write(6,27) n,iflag
      if( n .lt. 20 ) write(6,28) (x(j),j=1,n)
!                   done with this sample size
      end do  ! loop on i
      stop
      end program phsort
      SUBROUTINE HSORT(K,N)
!  HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N
!
!  ARGUMENTS
!  K   REAL VECTOR OF KEYS TO BE SORTED
!  N   NUMBER OF ITEMS TO BE SORTED
!
!  TO SORT A PARALLEL VECTOR OF RECORDS, USE HKSORT
!
!  J F MONAHAN (DEC, 1999) FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(N), INTENT(IN OUT) :: K
      REAL KK
      INTEGER I,L,NCUR
!
!                   DO NOTHING IF THERE'S NOTHING TO DO
      IF( N .LE. 1 ) RETURN
!                   INITIALIZE TO BUILDHEAP PART (LOOP ON L)
      L = N/2 + 1
      NCUR = N
      DO I = L,1,-1
      CALL HEAPIFY(I)
      END DO  ! LOOP ON I
      DO I = 2,N
!                   SWITCH CURRENT LARGEST WITH BOTTOM
      KK = K(1)
      K(1) = K(NCUR)
      K(NCUR) = KK
!                   REHEAP WITH ONE SHORTER
      NCUR = NCUR - 1
      CALL HEAPIFY(1)
      END DO  ! LOOP ON NCUR
      RETURN
      CONTAINS

      SUBROUTINE HEAPIFY(II)
      INTEGER, INTENT(IN) :: II
      INTEGER I,J
      I = II
      DO
      J = 2*I                 ! IS IT A LEAF OR ARE THERE SONS?
      IF( J > NCUR ) EXIT     ! A LEAF
      IF( J < NCUR ) THEN     ! ANOTHER SON OF I
                     IF( K(J+1) > K(J) ) J = J+1  ! LARGER SON IS K
                     END IF   ! A LEAF
      IF( K(J) > K(I) ) THEN  ! EXCHANGE
                        KK = K(J)
                        K(J) = K(I)
                        K(I) = KK
                        I = J
                        ELSE  ! EXIT -- HEAP PROPERTY
                        EXIT
                        END IF
      END DO ! WHILE NOT A LEAF
      END SUBROUTINE HEAPIFY

      END SUBROUTINE HSORT
      REAL FUNCTION RAN(IDUM)
!  PORTABLE IMPLEMENTATION OF UNIFORM PSEUDORANDOM NUMBER GENERATOR
!    LEWIS, GOODMAN, & MILLER MULTIPLICATIVE GENERATOR
!     X(N+1) = MOD( 16807*X(N), 2**31-1 )
!
!  P. BRANTLEY, B.L. FOX, L. SCHRAGE (1983) A GUIDE TO SIMULATION
!   SPRINGER-VERLAG, NEW YORK.  PP. 200-202.
!  UPDATED VERSION OF
!  LINUS SCHRAGE (1979)'A MORE PORTABLE FORTRAN RANDOM NUMBER GENERATOR' 
!   ACM TRANSACTIONS ON MATHEMATICAL SOFTWARE, VOLUME 5, PP. 132-138.
!
!  ARGUMENT 
!   IDUM    INTEGER  FIRST CALL SETS SEED, IGNORED IN SUBSEQUENT CALLS
!
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: IDUM
      REAL, PARAMETER :: RP = 4.656612875E-10     ! 1/P
      INTEGER, PARAMETER :: A = 16807             ! MULTIPLIER
      INTEGER, PARAMETER :: B = 127773            ! B = P / A
      INTEGER, PARAMETER :: C = 2836              ! C = P MOD A
      INTEGER, PARAMETER :: P = 2147483647        ! MODULUS 2**31 - 1
      INTEGER, SAVE :: IX = 0
      INTEGER K1
!
!                   IF NOT FIRST CALL, THEN SKIP SETTING SEED
      IF( IX .EQ. 0) IX = IDUM
!                   WRITE NUMBER AS ALPHA*2**16 + BETA
      K1 = IX / B
      IX = A*( IX - K1*B) - K1*C
!                   ABOVE DOES A*IX MOD B -K1*C
      IF( IX .LT. 0 ) IX = IX + P 
!                   RP BELOW IS RECIPROCAL OF P
      RAN = REAL(IX)*RP
      RETURN
      END FUNCTION RAN
! *** end of filename hsort.f95         *********************