1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
|
c subroutine quick (dataum, n, index)
c
c
c A quick-sorting algorithm suggested by the discussion on pages 114-119
c of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND SEARCHING, by
c D.E. Knuth, which was referenced in Don Wells' subroutine QUIK. This
c is my own attempt at encoding a quicksort-- PBS.
c
c Arguments
c
c datum (input / output) is a vector of dimension n containing randomly
c ordered real data upon input. Upon output the elements of
c dataum will be in order of increasing value.
c
c
c index (output) is an integer vector of dimension n. Upon return to
c the calling program the i-th element of index will tell where
c the i-th element of the sorted vector datum had been before
c datum was sorted.
c
c
c
c parameter (maxstack = 28)
c
c Parameters
c
c maxstack is the maximum number of entries the stack can contain.
c A limiting stack length of 28 restricts this quicksort
subroutine quick (datum, n, index, ier)
c
implicit none
integer n, index(n), ier
real datum(n)
c
real dkey
integer stklo(28), stkhi(28), i, lo, hi, nstak, limlo, limhi, ikey
c
c Initialize error code
ier = 0
c
c Initialize index.
c
do 50 i = 1, n
50 index(i) = i
c
c Initialize the pointers.
c
nstak = 0
limlo = 1
limhi = n
c
100 dkey = datum(limlo)
ikey = index(limlo)
c
c Compare all elements in the sub-vector between limlo and limhi with
c the current key datum.
c
lo = limlo
hi = limhi
101 continue
c
if (lo .eq. hi) go to 200
c
if (datum(hi) .le. dkey) go to 109
hi = hi - 1
c
c The pointer hi is to be left pointing at a datum smaller than the
c key, which is intended to be overwritten.
c
go to 101
c
109 datum(lo) = datum(hi)
index(lo) = index(hi)
lo = lo + 1
110 continue
c
if (lo .eq. hi) go to 200
c
if (datum(lo) .ge. dkey) go to 119
c
lo = lo + 1
go to 110
c
119 datum(hi) = datum(lo)
index(hi) = index(lo)
hi = hi - 1
c
c The pointer LO is to be left pointing at a datum LARGER than the
c key, which is intended to be overwritten.
c
go to 101
c
200 continue
c
c lo and hi are equal, and point at a value which is intended to
c be overwritten. Since all values below this point are less than
c the key and all values above this point are greater than the key,
c this is where we stick the key back into the vector.
c
datum(lo) = dkey
index(lo) = ikey
c
c At this point in the subroutine, all data between limlo and LO-1,
c inclusive, are less than datum(LO), and all data between LO+1 and
c limhi are larger than datum(LO).
c
c If both subarrays contain no more than one element, then take the most
c recent interval from the stack (if the stack is empty, we're done).
c If the larger of the two subarrays contains more than one element, and
c if the shorter subarray contains one or no elements, then forget the
c shorter one and reduce the other subarray. If the shorter subarray
c contains two or more elements, then place the larger subarray on the
c stack and process the subarray.
c
if ((limhi - lo) .gt. (lo - limlo)) go to 300
c
c Case 1: the lower subarray is longer. If it contains one or no
c elements then take the most recent interval from the stack and go
c back and operate on it.
c
if ((lo - limlo) .le. 1) go to 400
c
c If the upper (shorter) subinterval contains one or no elements, then
c process the lower (longer) one, but if the upper subinterval contains
c more than one element, then place the lower (longer) subinterval on
c the stack and process the upper one.
c
if ((limhi - lo) .ge. 2) go to 250
c
c Case 1a: the upper (shorter) subinterval contains no or one elements,
c so we go back and operate on the lower (longer) subinterval.
c
limhi = lo - 1
go to 100
c
250 continue
c
c Case 1b: the upper (shorter) subinterval contains at least two
c elements, so we place the lower (longer) subinterval on the stack and
c then go back and operate on the upper subinterval.
c
nstak = nstak + 1
if (nstak .gt. 28) then
ier = -1
return
endif
stklo(nstak) = limlo
stkhi(nstak) = lo - 1
limlo = lo + 1
go to 100
c
300 continue
c
c Case 2: the upper subarray is longer. If it contains one or no
c elements then take the most recent interval from the stack and
c operate on it.
c
if ((limhi - lo) .le. 1) go to 400
c
c If the lower (shorter) subinterval contains one or no elements, then
c process the upper (longer) one, but if the lower subinterval contains
c more than one element, then place the upper (longer) subinterval on
c the stack and process the lower one.
c
if ((lo - limlo) .ge. 2) go to 350
c
c Case 2a: the lower (shorter) subinterval contains no or one elements,
c so we go back and operate on the upper (longer) subinterval.
c
limlo = lo + 1
go to 100
c
350 continue
c
c Case 2b: the lower (shorter) subinterval contains at least two
c elements, so we place the upper (longer) subinterval on the stack and
c then go back and operate on the lower subinterval.
c
nstak = nstak + 1
if (nstak .gt. 28) then
ier = -1
return
endif
stklo(nstak) = lo + 1
stkhi(nstak) = limhi
limhi = lo - 1
go to 100
c
400 continue
c
c Take the most recent interval from the stack. If the stack happens
c to be empty, we are done.
c
if (nstak .le. 0) return
limlo = stklo(nstak)
limhi = stkhi(nstak)
nstak = nstak - 1
go to 100
c
end
|
