SUBROUTINE SORTRX(N,DATA,INDEX) C=================================================================== C C SORTRX -- SORT, Real input, indeX output C C C Input: N INTEGER C DATA REAL C C Output: INDEX INTEGER (DIMENSION N) C C This routine performs an in-memory sort of the first N elements of C array DATA, returning into array INDEX the indices of elements of C DATA arranged in ascending order. Thus, C C DATA(INDEX(1)) will be the smallest number in array DATA; C DATA(INDEX(N)) will be the largest number in DATA. C C The original data is not physically rearranged. The original order C of equal input values is not necessarily preserved. C C=================================================================== C C SORTRX uses a hybrid QuickSort algorithm, based on several C suggestions in Knuth, Volume 3, Section 5.2.2. In particular, the C "pivot key" [my term] for dividing each subsequence is chosen to be C the median of the first, last, and middle values of the subsequence; C and the QuickSort is cut off when a subsequence has 9 or fewer C elements, and a straight insertion sort of the entire array is done C at the end. The result is comparable to a pure insertion sort for C very short arrays, and very fast for very large arrays (of order 12 C micro-sec/element on the 3081K for arrays of 10K elements). It is C also not subject to the poor performance of the pure QuickSort on C partially ordered data. C C Created: 15 Jul 1986 Len Moss C C=================================================================== INTEGER N,INDEX(N) REAL DATA(N) INTEGER LSTK(31),RSTK(31),ISTK INTEGER L,R,I,J,P,INDEXP,INDEXT REAL DATAP C QuickSort Cutoff C C Quit QuickSort-ing when a subsequence contains M or fewer C elements and finish off at end with straight insertion sort. C According to Knuth, V.3, the optimum value of M is around 9. INTEGER M PARAMETER (M=9) C=================================================================== C C Make initial guess for INDEX DO 50 I=1,N INDEX(I)=I 50 CONTINUE C If array is short, skip QuickSort and go directly to C the straight insertion sort. IF (N.LE.M) GOTO 900 C=================================================================== C C QuickSort C C The "Qn:"s correspond roughly to steps in Algorithm Q, C Knuth, V.3, PP.116-117, modified to select the median C of the first, last, and middle elements as the "pivot C key" (in Knuth's notation, "K"). Also modified to leave C data in place and produce an INDEX array. To simplify C comments, let DATA[I]=DATA(INDEX(I)). C Q1: Initialize ISTK=0 L=1 R=N 200 CONTINUE C Q2: Sort the subsequence DATA[L]..DATA[R]. C C At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L, C r > R, and L <= m <= R. (First time through, there is no C DATA for l < L or r > R.) I=L J=R C Q2.5: Select pivot key C C Let the pivot, P, be the midpoint of this subsequence, C P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R) C so the corresponding DATA values are in increasing order. C The pivot key, DATAP, is then DATA[P]. P=(L+R)/2 INDEXP=INDEX(P) DATAP=DATA(INDEXP) IF (DATA(INDEX(L)) .GT. DATAP) THEN INDEX(P)=INDEX(L) INDEX(L)=INDEXP INDEXP=INDEX(P) DATAP=DATA(INDEXP) ENDIF IF (DATAP .GT. DATA(INDEX(R))) THEN IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN INDEX(P)=INDEX(L) INDEX(L)=INDEX(R) ELSE INDEX(P)=INDEX(R) ENDIF INDEX(R)=INDEXP INDEXP=INDEX(P) DATAP=DATA(INDEXP) ENDIF C Now we swap values between the right and left sides and/or C move DATAP until all smaller values are on the left and all C larger values are on the right. Neither the left or right C side will be internally ordered yet; however, DATAP will be C 300 CONTINUE C Q3: Search for datum on left >= DATAP C C At this point, DATA[L] <= DATAP. We can therefore start scanning C up from L, looking for a value >= DATAP (this scan is guaranteed C to terminate since we initially placed DATAP near the middle of C the subsequence). I=I+1 IF (DATA(INDEX(I)).LT.DATAP) GOTO 300 400 CONTINUE C Q4: Search for datum on right <= DATAP C C At this point, DATA[R] >= DATAP. We can therefore start scanning C down from R, looking for a value <= DATAP (this scan is guaranteed C to terminate since we initially placed DATAP near the middle of C the subsequence). J=J-1 IF (DATA(INDEX(J)).GT.DATAP) GOTO 400 C Q5: Have the two scans collided? IF (I.LT.J) THEN C Q6: No, interchange DATA[I] <--> DATA[J] and continue INDEXT=INDEX(I) INDEX(I)=INDEX(J) INDEX(J)=INDEXT GOTO 300 ELSE C Q7: Yes, select next subsequence to sort C C At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r], C for all L <= l < I and J < r <= R. If both subsequences are C more than M elements long, push the longer one on the stack and C go back to QuickSort the shorter; if only one is more than M C elements long, go back and QuickSort it; otherwise, pop a C subsequence off the stack and QuickSort it. IF (R-J .GE. I-L .AND. I-L .GT. M) THEN ISTK=ISTK+1 LSTK(ISTK)=J+1 RSTK(ISTK)=R R=I-1 ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN ISTK=ISTK+1 LSTK(ISTK)=L RSTK(ISTK)=I-1 L=J+1 ELSE IF (R-J .GT. M) THEN L=J+1 ELSE IF (I-L .GT. M) THEN R=I-1 ELSE C Q8: Pop the stack, or terminate QuickSort if empty IF (ISTK.LT.1) GOTO 900 L=LSTK(ISTK) R=RSTK(ISTK) ISTK=ISTK-1 ENDIF GOTO 200 ENDIF 900 CONTINUE C=================================================================== C C Q9: Straight Insertion sort DO 950 I=2,N IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN INDEXP=INDEX(I) DATAP=DATA(INDEXP) P=I-1 920 CONTINUE INDEX(P+1) = INDEX(P) P=P-1 IF (P.GT.0) THEN IF (DATA(INDEX(P)).GT.DATAP) GOTO 920 ENDIF INDEX(P+1) = INDEXP ENDIF 950 CONTINUE C=================================================================== C C All done END