!     Last change:  DOS   2 Aug 2000    9:14 pm
!                                       *** copyright 2000 ***
! *** filename dks2s.f95                *** John F. Monahan **
!                                       **********************
      program pdks2s
!   test of dks2s -- computes Kolmogorov-Smirnov two sample statistic
!
      implicit none
      real, dimension(100) :: x,y,z
      real s,dp,dn,dpp,dnn
      real ran
      integer, dimension(100) :: list
      integer  m,n,mpn,i,il,si,kprob
  20  FORMAT(//' Textbook Test Problem with Sample Sizes',2i4)
  21  FORMAT(//' Generated Test Problem',i3,' with Sample Sizes',2i4)
  22  FORMAT(/' K-S statistics computed with DKS2S',2f8.4)
  23  FORMAT(/' K-S statistics computed in textbook fashion',2f8.4)
  26  FORMAT(2F7.2)
  27  FORMAT(12f6.2)
!
!                   open output file
      open( unit=6, file='dks2s.out' )
!                   open data file
      open( unit=5, file='dks2s.dat' )
!                   initialize uniform generator
      x(1) = ran(5151917)
!                   first a case out of Hollander & Wolfe
!                     data attributed to Delse & Feather
      m = 10
      n = 10
      write(6,20) m,n
      do i = 1,10
      read(5,26) x(i),y(i)
      write(6,26) x(i),y(i)
      end do  ! loop on i
!                   use dks2s to compute Kolmogorov-Smirnov
!                     two sample statistic (both sides)
      call hsort(x,m)
      call hsort(y,n)
      call dks2s(x,m,y,n,dp,dn)
      write(6,22) dp,dn
!                   let's compute it another way
      mpn = m + n
      do i = 1,m
      z(i) = x(i)
      list(i) = 1
      end do  ! loop on i
      do i = 1,n
      z(m+i) = y(i)
      list(m+i) = 0
      end do  ! loop on i
!                   sort the joint vector
!                     along with list
      call hksort(z,list,mpn)
!                   compute Kolmogorov-Smirnov two sample statistic
!                     the textbook way
      dpp = 0.
      dnn = 0.
      s = 0.
      i = 1
      il = 0
      si = 0
      do while( i + il .le. mpn )
      do while ( z(i) .eq. z(i+il) )    ! collect ties
      si = si + list(i+il) 
      il = il + 1
      if( i + il .GT. mpn ) exit
      end do  ! while( z(i) .eq. z(i+il)
      s = s + float(si)/float(m) - float(il-si)/float(n)
      dpp = max( dpp, s )
      dnn = min( dnn, s )
      i = i + il
      il = 0
      si = 0
      end do  ! while( i + il .le. mpn )
      write(6,23) dpp,dnn
!
!                        *** now some randomly generated problems ***
!
      do kprob = 1,20
      m = 1 + int( 10.*ran(kprob) )
      n = 1 + int( 40.*ran(kprob) )
      write(6,21) kprob,m,n
!                   generate data so that it has some ties
      do i = 1,m
      x(i) = int( 20.*(ran(i)-.3) ) / 10.
      z(i) = x(i)
      list(i) = 1
      end do  ! loop on i
      do i = 1,n
      y(i) = int( 30.*(ran(i)-.4) ) / 10.
      z(i+m) = y(i)
      list(i+m) = 0
      end do  ! loop on i
      write(6,27) (x(i),i=1,m)
      write(6,27) (y(i),i=1,n)
!                   use dks2s to compute Kolmogorov-Smirnov
!                     two sample statistic (both sides)
      call hsort(x,m)
      call hsort(y,n)
      call dks2s(x,m,y,n,dp,dn)
      write(6,22) dp,dn
!                   let's compute it another way
      mpn = m + n
!                   sort the joint vector
!                     along with list
      call hksort(z,list,mpn)
!                   compute Kolmogorov-Smirnov two sample statistic
!                     the textbook way
      dpp = 0.
      dnn = 0.
      s = 0.
      i = 1
      il = 0
      si = 0
      do while( i + il .le. mpn )
      do while ( z(i) .eq. z(i+il) )    ! collect ties
      si = si + list(i+il) 
      il = il + 1
      if( i + il .GT. mpn ) exit
      end do  ! while( z(i) .eq. z(i+il)
      s = s + real(si)/real(m) - real(il-si)/real(n)
      dpp = max( dpp, s )
      dnn = min( dnn, s )
      i = i + il
      il = 0
      si = 0
      end do  ! while( i + il .le. mpn )
      write(6,23) dpp,dnn
      end do  ! loop on kprob
      stop
      end program pdks2s
      subroutine dks2s(x,m,y,n,dp,dn)
!  computes Kolmogorov-Smirnov two sample test statistics
!
!  *** requires both x and y to be already sorted ***
!
!  x   real vector of m observations
!  m   integer    number of x's
!  y   real vector of n observations
!  n   integer    number of y's
!  dp  real       distance in the positive direction
!  dn  real       distance in the negative direction
!
! J F Monahan (June 1998) Dept of Statistics, NC State Univ, Raleigh, USA
!    recoded February, April 2000 for Fortran 95
!
      implicit none
      integer, intent(in) :: m,n
      real, dimension(m), intent(in) :: x
      real, dimension(n), intent(in) :: y
      real, intent(out) :: dp,dn
      integer j,jl,k,kl,s,sp,sn
      s = 0
      sp = 0
      sn = 0
!                   complete initializations
      j = 1
      k = 1
      jl = 0
      kl = 0
!                   start/restart here
      do  ! unrestricted do
!                   how many x's are the same?
      if( jl .eq. 0 ) then
      do while( x(j) .eq. x(j+jl) )
      jl = jl + 1
      if( j + jl .gt. m ) exit
      end do  ! while tied
                      end if ! ( jl .eq. 0 )
!                   have jl x's that are equal
!                   how many y's are the same?
      if( kl .eq. 0 ) then
      do while( y(k) .eq. y(k+kl) )
      kl = kl + 1
      if( k + kl .gt. n ) exit
      end do  ! while tied
                       end if ! ( kl .eq. 0 )
!                   have kl y's that are equal
!                   now which are bigger?
      if( x(j) .eq. y(k) ) then
!                   here x(j) = y(k)
      s = s + jl*n - kl*m
      sp = max( sp, s )
      sn = min( sn, s )
!                   hit either boundary and we're done
      if( j + jl .gt. m ) exit
      if( k + kl .gt. n ) exit
      j = j + jl
      jl = 0
      k = k + kl
      kl = 0        ! end of x and y are equal
                           else         ! unequal
      if( x(j) .lt. y(k) ) then
!                    x is smaller
      s = s + jl*n
      sp = max( sp, s )
!                   hit either boundary and we're done
      if( j + jl .gt. m ) exit
      j = j + jl
      jl = 0
                           else
!                   y is smaller
      s = s - kl*m
      sn = min( sn, s )
!                   hit either boundary and we're done
      if( k + kl .gt. n ) exit
      k = k + kl
      kl = 0
                           end if ! x is smaller
                           end if ! x and y are equal
      end do ! unrestricted do
      dp = real(sp)/real(m*n)
      dn = real(sn)/real(m*n)
      return
      end subroutine dks2s
      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
      SUBROUTINE HKSORT(K,M,N)
!  HEAPSORT ALGORITHM FOR SORTING ON VECTOR OF KEYS K OF LENGTH N
!
!  ARGUMENTS
!  K   REAL VECTOR OF KEYS TO BE SORTED
!  M   INTEGER VECTOR TO BE MOVED IN PARALLEL TO K
!  N   NUMBER OF ITEMS TO BE SORTED
!
!  TO SORT A PARALLEL VECTOR OF RECORDS, USE THIS ONE
!
!  J F MONAHAN (DEC, 1999) FORTRAN 95
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: N
      REAL, DIMENSION(N), INTENT(IN OUT) :: K
      INTEGER, DIMENSION(N), INTENT(IN OUT) :: M
      REAL KK
      INTEGER I,L,NCUR,MM
!
!                   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
!                   ALSO SWITCH PARALLEL VECTOR
      MM = M(1)
      M(1) = M(NCUR)
      M(NCUR) = MM
!                   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
!                   ALSO SWITCH PARALLEL VECTOR
      MM = M(J)
      M(J) = M(I)
      M(I) = MM
                        I = J
                        ELSE  ! EXIT -- HEAP PROPERTY
                        EXIT
                        END IF
      END DO ! WHILE NOT A LEAF
      END SUBROUTINE HEAPIFY

      END SUBROUTINE HKSORT
! *** end of filename dks2s.f95         *********************