! Last change: DOS 3 Aug 2000 6:12 pm ! *** copyright 2000 *** ! *** filename partit.f95 *** John F. Monahan ** ! ********************** ! TEST OF PARTITION ALGORITHM PARTIT program ppartit implicit none real, dimension(3600) :: x real xsave integer i,j,k,n,ns,nsin,nt,iflag real ran ! 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) 24 FORMAT(' NSIN=',I4,F12.6,' NS,NT ',2I4) 26 FORMAT(2X,8F9.6) 28 FORMAT(2X,12F6.2) ! open output file open( unit=6, file='partit.out' ) ! initialize uniform generator x(1) = ran(5151917) ! let's check out these twenty cases do i = 1,20 n = i if( i .gt. 10 ) n = i*(i-8) write(6,20) n ! first all positive do j = 1,n x(j) = 5.*ran(j) end do ! loop on j if( n .lt. 20 ) write(6,22) (x(j),j=1,n) ! select random element to be used to partition ns = 1 + int( ran(1)*real(n) ) xsave = x(ns) nsin = ns ! call partitioning algorithm and write out results call partit(x,1,n,ns,nt) write(6,24) nsin,xsave,ns,nt ! test iflag = 0 do j = 1,ns if( x(j) .ge. xsave ) iflag = 1 end do ! loop on j do j = ns+1,ns+nt if( x(j) .ne. xsave ) iflag = 1 end do ! loop on j do j = ns+nt+1,n if( x(j) .le. xsave ) 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 if( n .lt. 20 ) write(6,22) (x(j),j=1,n) ! select random element to be used to partition ns = 1 + int( ran(1)*real(n) ) xsave = x(ns) nsin = ns ! call partitioning algorithm and write out results call partit(x,1,n,ns,nt) WRITE(6,24) nsin,xsave,ns,nt ! test iflag = 0 do j = 1,ns if( x(j) .ge. xsave ) iflag = 1 end do ! loop on j do j = ns+1,ns+nt if( x(j) .ne. xsave ) iflag = 1 end do ! loop on j do j = ns+nt+1,n if( x(j) .le. xsave ) 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 if( n .lt. 20 ) write(6,26) (x(j),j=1,n) ! select random element to be used to partition ns = 1 + int( ran(1)*real(n) ) xsave = x(ns) nsin = ns ! call partitioning algorithm and write out results call partit(x,1,n,ns,nt) WRITE(6,24) nsin,xsave,ns,nt ! test iflag = 0 do j = 1,ns if( x(j) .ge. xsave ) iflag = 1 end do ! loop on j do j = ns+1,ns+nt if( x(j) .ne. xsave ) iflag = 1 end do ! loop on j do j = ns+nt+1,n if( x(j) .le. xsave ) 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 = 8.*(ran(j)-.3) x(j) = float(k)/10. end do ! loop on j if( n .lt. 20 ) write(6,28) (x(j),j=1,n) ! select random element to be used to partition ns = 1 + int( ran(1)*real(n) ) xsave = x(ns) nsin = ns ! call partitioning algorithm and write out results call partit(x,1,n,ns,nt) write(6,24) nsin,xsave,ns,nt ! test iflag = 0 do j = 1,ns if( x(j) .ge. xsave ) iflag = 1 end do ! loop on j do j = ns+1,ns+nt if( x(j) .ne. xsave ) iflag = 1 end do ! loop on j do j = ns+nt+1,n if( x(j) .le. xsave ) 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 ppartit SUBROUTINE PARTIT(X,LEFT,RGHT,NS,NT) ! ALGORITHM FOR PARTITIONING LIST OF KEYS K FROM X(LEFT) TO X(RGHT) ! INTO THREE GROUPS: THOSE < X*, THOSE = X*, THOSE > X* ! ! ARGUMENTS ! X IN REAL LIST OF KEYS TO BE PARTITIONED ! OUT PARTITIONED LIST ! X(LEFT) ... X(LEFT+NS-1) < X* ! X(LEFT+NS) ... X(LEFT+NS+NT-1) = X* ! X(LEFT+NS+NT) ... X(RGHT) > X* ! ! LEFT IN INTEGER LEFT ENDPOINT OF LIST ! RGHT IN INTEGER RIGHT ENDPOINT OF LIST ! NS IN INTEGER POINTS TO PARTITION ELEMENT X* = X(NS) (IN) ! OUT NUMBER IN LIST < X* ! NT OUT INTEGER NUMBER IN LIST = X* ! IMPLICIT NONE INTEGER, INTENT(IN) :: LEFT,RGHT INTEGER, INTENT(IN OUT) :: NS INTEGER, INTENT(OUT) :: NT REAL, DIMENSION(RGHT) :: X REAL XCUT,XX INTEGER I,J ! USE X(NS) TO SPLIT X XCUT = X(NS) ! DO NOTHING IF THERE'S NOTHING TO DO NS = 0 NT = 1 IF( RGHT - LEFT .LE. 0 ) RETURN ! FIRST PARITITION TO TWO GROUPS: LT AND GE ! ! INITIALIZE POINTERS I (LEFT) AND J (RIGHT) I = LEFT J = RGHT DO WHILE( I .LE. J ) DO WHILE( X(I) .LT. XCUT ) I = I + 1 IF( I .GT. J ) EXIT ! LAST ONE LT IS I-1 END DO ! WHILE ( X(I) .LT. XCUT ) IF( I .GT. J ) EXIT ! FINISH ! HAVE X(I) OUT OF PLACE DO WHILE ( X(J) .GE. XCUT ) J = J - 1 IF( I .GE. J ) EXIT END DO ! WHILE ( X(J) .GE. XCUT ) IF( I .GE. J ) EXIT ! HAVE X(J) OUT OF PLACE ! EXCHANGE I AND J XX = X(I) X(I) = X(J) X(J) = XX I = I + 1 J = J - 1 END DO ! FIRST MAIN WHILE ! HALFWAY DONE -- X(I-1) IS END OF SMALLER NS = I - LEFT ! NOW SEPARATE THE TIED FROM THE LARGER ! ! ! INITIALIZE POINTERS I (LEFT) AND J (RIGHT) J = RGHT DO WHILE( I .LE. J ) DO WHILE( X(I) .EQ. XCUT ) I = I + 1 IF( I .GT. J ) EXIT ! LAST ONE EQ IS I-1 END DO ! WHILE ( X(I) .EQ. XCUT ) IF( I .GT. J ) EXIT ! FINISH ! HAVE X(I) OUT OF PLACE DO WHILE( X(J) .GT. XCUT ) J = J - 1 IF( I .GE. J ) EXIT END DO ! WHILE ( X(J) .GT. XCUT ) IF( I .GE. J ) EXIT ! HAVE X(J) OUT OF PLACE ! EXCHANGE I AND J XX = X(I) X(I) = X(J) X(J) = XX I = I + 1 J = J - 1 END DO ! SECOND MAIN WHILE NT = I - LEFT - NS RETURN END SUBROUTINE PARTIT 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 partit.f95 *********************