quick_sort2.f 6.8 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232
  1. C From Leonard J. Moss of SLAC:
  2. C Here's a hybrid QuickSort I wrote a number of years ago. It's
  3. C based on suggestions in Knuth, Volume 3, and performs much better
  4. C than a pure QuickSort on short or partially ordered input arrays.
  5. SUBROUTINE SORTRX(N,DATA,INDEX)
  6. C===================================================================
  7. C
  8. C SORTRX -- SORT, Real input, indeX output
  9. C
  10. C
  11. C Input: N INTEGER
  12. C DATA REAL
  13. C
  14. C Output: INDEX INTEGER (DIMENSION N)
  15. C
  16. C This routine performs an in-memory sort of the first N elements of
  17. C array DATA, returning into array INDEX the indices of elements of
  18. C DATA arranged in ascending order. Thus,
  19. C
  20. C DATA(INDEX(1)) will be the smallest number in array DATA;
  21. C DATA(INDEX(N)) will be the largest number in DATA.
  22. C
  23. C The original data is not physically rearranged. The original order
  24. C of equal input values is not necessarily preserved.
  25. C
  26. C===================================================================
  27. C
  28. C SORTRX uses a hybrid QuickSort algorithm, based on several
  29. C suggestions in Knuth, Volume 3, Section 5.2.2. In particular, the
  30. C "pivot key" [my term] for dividing each subsequence is chosen to be
  31. C the median of the first, last, and middle values of the subsequence;
  32. C and the QuickSort is cut off when a subsequence has 9 or fewer
  33. C elements, and a straight insertion sort of the entire array is done
  34. C at the end. The result is comparable to a pure insertion sort for
  35. C very short arrays, and very fast for very large arrays (of order 12
  36. C micro-sec/element on the 3081K for arrays of 10K elements). It is
  37. C also not subject to the poor performance of the pure QuickSort on
  38. C partially ordered data.
  39. C
  40. C Created: 15 Jul 1986 Len Moss
  41. C
  42. C===================================================================
  43. INTEGER N,INDEX(N)
  44. REAL*8 DATA(N)
  45. INTEGER LSTK(31),RSTK(31),ISTK
  46. INTEGER L,R,I,J,P,INDEXP,INDEXT
  47. REAL*8 DATAP
  48. C QuickSort Cutoff
  49. C
  50. C Quit QuickSort-ing when a subsequence contains M or fewer
  51. C elements and finish off at end with straight insertion sort.
  52. C According to Knuth, V.3, the optimum value of M is around 9.
  53. INTEGER M
  54. PARAMETER (M=9)
  55. C===================================================================
  56. C
  57. C Make initial guess for INDEX
  58. DO 50 I=1,N
  59. INDEX(I)=I
  60. 50 CONTINUE
  61. C If array is short, skip QuickSort and go directly to
  62. C the straight insertion sort.
  63. IF (N.LE.M) GOTO 900
  64. C===================================================================
  65. C
  66. C QuickSort
  67. C
  68. C The "Qn:"s correspond roughly to steps in Algorithm Q,
  69. C Knuth, V.3, PP.116-117, modified to select the median
  70. C of the first, last, and middle elements as the "pivot
  71. C key" (in Knuth's notation, "K"). Also modified to leave
  72. C data in place and produce an INDEX array. To simplify
  73. C comments, let DATA[I]=DATA(INDEX(I)).
  74. C Q1: Initialize
  75. ISTK=0
  76. L=1
  77. R=N
  78. 200 CONTINUE
  79. C Q2: Sort the subsequence DATA[L]..DATA[R].
  80. C
  81. C At this point, DATA[l] <= DATA[m] <= DATA[r] for all l < L,
  82. C r > R, and L <= m <= R. (First time through, there is no
  83. C DATA for l < L or r > R.)
  84. I=L
  85. J=R
  86. C Q2.5: Select pivot key
  87. C
  88. C Let the pivot, P, be the midpoint of this subsequence,
  89. C P=(L+R)/2; then rearrange INDEX(L), INDEX(P), and INDEX(R)
  90. C so the corresponding DATA values are in increasing order.
  91. C The pivot key, DATAP, is then DATA[P].
  92. P=(L+R)/2
  93. INDEXP=INDEX(P)
  94. DATAP=DATA(INDEXP)
  95. IF (DATA(INDEX(L)) .GT. DATAP) THEN
  96. INDEX(P)=INDEX(L)
  97. INDEX(L)=INDEXP
  98. INDEXP=INDEX(P)
  99. DATAP=DATA(INDEXP)
  100. ENDIF
  101. IF (DATAP .GT. DATA(INDEX(R))) THEN
  102. IF (DATA(INDEX(L)) .GT. DATA(INDEX(R))) THEN
  103. INDEX(P)=INDEX(L)
  104. INDEX(L)=INDEX(R)
  105. ELSE
  106. INDEX(P)=INDEX(R)
  107. ENDIF
  108. INDEX(R)=INDEXP
  109. INDEXP=INDEX(P)
  110. DATAP=DATA(INDEXP)
  111. ENDIF
  112. C Now we swap values between the right and left sides and/or
  113. C move DATAP until all smaller values are on the left and all
  114. C larger values are on the right. Neither the left or right
  115. C side will be internally ordered yet; however, DATAP will be
  116. C in its final position.
  117. 300 CONTINUE
  118. C Q3: Search for datum on left >= DATAP
  119. C
  120. C At this point, DATA[L] <= DATAP. We can therefore start scanning
  121. C up from L, looking for a value >= DATAP (this scan is guaranteed
  122. C to terminate since we initially placed DATAP near the middle of
  123. C the subsequence).
  124. I=I+1
  125. IF (DATA(INDEX(I)).LT.DATAP) GOTO 300
  126. 400 CONTINUE
  127. C Q4: Search for datum on right <= DATAP
  128. C
  129. C At this point, DATA[R] >= DATAP. We can therefore start scanning
  130. C down from R, looking for a value <= DATAP (this scan is guaranteed
  131. C to terminate since we initially placed DATAP near the middle of
  132. C the subsequence).
  133. J=J-1
  134. IF (DATA(INDEX(J)).GT.DATAP) GOTO 400
  135. C Q5: Have the two scans collided?
  136. IF (I.LT.J) THEN
  137. C Q6: No, interchange DATA[I] <--> DATA[J] and continue
  138. INDEXT=INDEX(I)
  139. INDEX(I)=INDEX(J)
  140. INDEX(J)=INDEXT
  141. GOTO 300
  142. ELSE
  143. C Q7: Yes, select next subsequence to sort
  144. C
  145. C At this point, I >= J and DATA[l] <= DATA[I] == DATAP <= DATA[r],
  146. C for all L <= l < I and J < r <= R. If both subsequences are
  147. C more than M elements long, push the longer one on the stack and
  148. C go back to QuickSort the shorter; if only one is more than M
  149. C elements long, go back and QuickSort it; otherwise, pop a
  150. C subsequence off the stack and QuickSort it.
  151. IF (R-J .GE. I-L .AND. I-L .GT. M) THEN
  152. ISTK=ISTK+1
  153. LSTK(ISTK)=J+1
  154. RSTK(ISTK)=R
  155. R=I-1
  156. ELSE IF (I-L .GT. R-J .AND. R-J .GT. M) THEN
  157. ISTK=ISTK+1
  158. LSTK(ISTK)=L
  159. RSTK(ISTK)=I-1
  160. L=J+1
  161. ELSE IF (R-J .GT. M) THEN
  162. L=J+1
  163. ELSE IF (I-L .GT. M) THEN
  164. R=I-1
  165. ELSE
  166. C Q8: Pop the stack, or terminate QuickSort if empty
  167. IF (ISTK.LT.1) GOTO 900
  168. L=LSTK(ISTK)
  169. R=RSTK(ISTK)
  170. ISTK=ISTK-1
  171. ENDIF
  172. GOTO 200
  173. ENDIF
  174. 900 CONTINUE
  175. C===================================================================
  176. C
  177. C Q9: Straight Insertion sort
  178. DO 950 I=2,N
  179. IF (DATA(INDEX(I-1)) .GT. DATA(INDEX(I))) THEN
  180. INDEXP=INDEX(I)
  181. DATAP=DATA(INDEXP)
  182. P=I-1
  183. 920 CONTINUE
  184. INDEX(P+1) = INDEX(P)
  185. P=P-1
  186. IF (P.GT.0) THEN
  187. IF (DATA(INDEX(P)).GT.DATAP) GOTO 920
  188. ENDIF
  189. INDEX(P+1) = INDEXP
  190. ENDIF
  191. 950 CONTINUE
  192. C===================================================================
  193. C
  194. C All done
  195. END