git.iko.soy/bat
on master help sign in
Add Fortran (Fixed Form) syntax test file
commitb02120cf6679b50bb4e4b2dd9c8e72be5ed4372b authorMohamed Abdelnour <mohamed.k.abdelnour@gmail.com> dateSun 30 May 2021, 11:42 (5 years ago) parent702b5ca
3 files changed+671
compare with master
unifiedsplit
tests/syntax-tests/highlighted/Fortran (Fixed Form)/quicksort_real_F77.F+323 tests/syntax-tests/source/Fortran (Fixed Form)/LICENSE.md+25 tests/syntax-tests/source/Fortran (Fixed Form)/quicksort_real_F77.F+323
tests/syntax-tests/highlighted/Fortran (Fixed Form)/quicksort_real_F77.Fnew file+323
@@ -0,0 +1,323 @@
1+C Fortran 77 implementation of a quicksort algorithm for arrays with
2+C real entries.
3+C ----------
4+C June 2019 
5+C Jason Allen Anema, Ph.D.
6+C Division of Statistical Genomics
7+C Department of Genetics
8+C Washington University School of Medicine in St. Louis
9+C 
10+C This work is partially supported NIH grant AG023746
11+C ----------
12+C Insertion sort is used for short arrays, as quicksort is slower on
13+C these.
14+C 
15+C Hoare partition scheme is used (sweeping left and right), as it does
16+C three times fewer swaps on average that the Lamuto partition scheme.
17+C In conjunction with this, tripartite partition is performed
18+C concurrently (solving the "Dutch National Flag problem"). This avoids 
19+C horrible runtimes on highly repetitive arrays. For example, without 
20+C this, an array of random zeros and ones would have a runtime of
21+C O(N^2), but now has a runtime of O(N). The runtime for this algorthm
22+C on arrays with k highly repetitive entries is now O(kN).
23+C 
24+C For medium length (sub)arrays, pivots are choosen using
25+C Median-of-Three, and those three items are sorted. For longer (sub)arrays
26+C the pseudomedian of nine (Median of medians). This avoids O(N^2) runtime on
27+C nonrandom inputs such as increasing and decreasing sequences. 
28+C
29+C See Louis Bentley, Jon & McIlroy, Douglas. (1993). Engineering a Sort Function.
30+C Softw., Pract. Exper.. 23. 1249-1265. 10.1002/spe.4380231105 for details. 
31+C
32+C The ordering on elements of the array are defined by a comparison
33+C function,compar, that is a user-supplied INTEGER*2 function of the form
34+C compar(a,b) which returns:
35+C -1 if a precedes b
36+C +1 if b precedes a
37+C 0 is a and b are considered equivalent
38+C and thus defines a total ordering.
39+C 
40+C If one would like to use the standard order on integers, the
41+C compar function could be written in a file "compint.F" as:
42+C ----------------------------------------------------------------
43+C INTEGER*2 FUNCTION compint(a,b)
44+C INTEGER a, b
45+C if(a.lt.b)then
46+C compint = -1
47+C elseif(a.gt.b)then
48+C compint = +1
49+C else
50+C compint = 0
51+C endif
52+C END
53+C ----------------------------------------------------------------
54+C Then in your program, call quicksort with:
55+C call quicksort_real_F77(array, n, compint)
56+C
57+C The maximal length of an array in this implementation is (2^31-1),
58+C but can be changed to allow for length up to (2^63-1) by changing the
59+C data types of the relevant variables and constants. If you wish to 
60+C sort longer arrays, of length N, you'll need to customize variable 
61+C and constant types and set mstack to be at least (2*log_2(N)+2). 
62+C
63+C ----------------------------------------------------------------
64+C Copyright 2019 Jason Allen Anema
65+C 
66+C Permission is hereby granted, free of charge, to any person obtaining
67+C a copy of this software and associated documentation files (the "Software"),
68+C to deal in the Software without restriction, including without limitation the
69+C rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
70+C sell copies of the Software, and to permit persons to whom the Software is
71+C furnished to do so, subject to the following conditions:
72+C
73+C The above copyright notice and this permission notice shall be included
74+C in all copies or substantial portions of the Software.
75+C
76+C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
77+C OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
78+C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
79+C THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
80+C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
81+C FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
82+C IN THE SOFTWARE.
83+C -------------------------------------------------------------------
84+C
85+ SUBROUTINE quicksort_real_F77(array,n,compar)
86+ INTEGER n, maxins, maxmid, mstack
87+ REAL array(n)
88+ PARAMETER (maxins = 7, maxmid = 40, mstack = 128)
89+C maxins: maximal size of (sub)arrays to be sorted with
90+C insertion sort.
91+C maxmid: maximal size of (sub)arrays that will be quicksorted with
92+C Median-of-Three pivots.
93+C mstack: maximal size of required auxiliary storage (a stack), plus 2 
94+C extra spots, which tracks the starts and ends of yet unsorted 
95+C subarrays. mstack = 130 is large enough to handle arrays up to 
96+C length 2^63-1. This maximal size follows from
97+C processing smaller arrays first and pigeonhole principal.
98+C 
99+ INTEGER a, d, i, j, k, s, lo, mid, hi, tstack, bstack(mstack)
100+C a, d, i, j, k, s: indices
101+C lo, mid, and hi: their natural location in a (sub)array
102+C tstack: equal to twice the number of additional subarrays still 
103+C needing to be sorted
104+C bstack: stack of the endpoints of unsorted subarrays
105+ INTEGER pm1, pm2, pm3, pm4, pm5, pm6, pm7, pm8, pm9
106+C for pseudomedian of nine positions in (sub)arrays
107+ REAL piv, temp
108+C piv is to store the pivot's value
109+C 
110+ EXTERNAL compar
111+ INTEGER*2 compar
112+C compar is a user-supplied INTEGER*2 function of the form
113+C compar(a,b) which returns:
114+C -1 if a precedes b
115+C +1 if b precedes a
116+C 0 is a and b are considered equivalent
117+C and thus defines a total ordering. 
118+ tstack = 0
119+ lo = 1 
120+ hi = n
121+C
122+C Insertion sort subarrays of size maxins or less
123+ 1 if(hi-lo+1.le.maxins)then
124+ do 10, i = lo + 1, hi, 1
125+ temp = array(i)
126+ do 11 j = i - 1, lo, -1
127+ if(compar(array(j), temp).le.0)goto 2
128+ array(j+1)=array(j)
129+ 11 continue
130+ j = lo - 1
131+ 2 array(j+1) = temp
132+ 10 continue
133+ if(tstack.eq.0)return
134+C Pop the bstack, and start new partitioning
135+ hi = bstack(tstack)
136+ lo = bstack(tstack-1)
137+ tstack = tstack - 2
138+ else
139+C Use Median-of-Three as choice of pivot (median of lo, middle, hi)
140+C and reorder those elements appropriately when subarrays are medium
141+C length (between maxins and maxmid)
142+ mid = lo + (hi-lo)/2
143+ if(hi-lo.le.maxmid)then
144+ if(compar(array(mid), array(lo)).eq.-1)then
145+ temp = array(lo)
146+ array(lo) = array(mid)
147+ array(mid) = temp
148+ endif
149+ if(compar(array(hi), array(lo)).eq.-1)then
150+ temp = array(hi)
151+ array(hi) = array(lo)
152+ array(lo) = temp
153+ endif
154+ if(compar(array(hi), array(mid)).eq.-1)then
155+ temp = array(hi)
156+ array(hi) = array(mid)
157+ array(mid) = temp
158+ endif
159+C Use pseudomedian of nine (Median of medians) as choice of pivot when 
160+C subarrays are longer than maxmid. Note that doing it this way requires only 12
161+C comparisons for finding the pivot.
162+ elseif(hi-lo+1.gt.maxmid)then
163+ pm1 = lo
164+ pm5 = lo + (hi-lo)/2
165+ pm9 = hi
166+ pm3 = lo + (pm5-lo)/2
167+ pm7 = pm5 + (hi-pm5)/2
168+ pm2 = lo + (pm3-lo)/2
169+ pm4 = pm3 + (pm5-pm3)/2
170+ pm6 = pm5 + (pm7-pm5)/2
171+ pm8 = pm7 + (pm9-pm7)/2
172+C Median and sorting for pm1, pm2, pm3
173+ if(compar(array(pm2), array(pm1)).eq.-1)then
174+ temp = array(pm1)
175+ array(pm1) = array(pm2)
176+ array(pm2) = temp
177+ endif
178+ if(compar(array(pm3), array(pm1)).eq.-1)then
179+ temp = array(pm3)
180+ array(pm3) = array(pm1)
181+ array(pm1) = temp
182+ endif
183+ if(compar(array(pm3), array(pm2)).eq.-1)then
184+ temp = array(pm3)
185+ array(pm3) = array(pm2)
186+ array(pm2) = temp
187+ endif
188+C Median and sorting for pm4, pm5, pm6
189+ if(compar(array(pm5), array(pm4)).eq.-1)then
190+ temp = array(pm4)
191+ array(pm4) = array(pm5)
192+ array(pm5) = temp
193+ endif
194+ if(compar(array(pm6), array(pm4)).eq.-1)then
195+ temp = array(pm6)
196+ array(pm6) = array(pm4)
197+ array(pm4) = temp
198+ endif
199+ if(compar(array(pm6), array(pm5)).eq.-1)then
200+ temp = array(pm6)
201+ array(pm6) = array(pm5)
202+ array(pm5) = temp
203+ endif
204+C Median and sorting for pm7, pm8, pm9
205+ if(compar(array(pm8), array(pm7)).eq.-1)then
206+ temp = array(pm7)
207+ array(pm7) = array(pm8)
208+ array(pm8) = temp
209+ endif
210+ if(compar(array(pm9), array(pm7)).eq.-1)then
211+ temp = array(pm9)
212+ array(pm9) = array(pm7)
213+ array(pm7) = temp
214+ endif
215+ if(compar(array(pm9), array(pm8)).eq.-1)then
216+ temp = array(pm9)
217+ array(pm9) = array(pm8)
218+ array(pm8) = temp
219+ endif
220+C Median of the medians (which are now pm2, pm5, pm8)
221+ if(compar(array(pm5), array(pm2)).eq.-1)then
222+ temp = array(pm2)
223+ array(pm2) = array(pm5)
224+ array(pm5) = temp
225+ endif
226+ if(compar(array(pm8), array(pm2)).eq.-1)then
227+ temp = array(pm8)
228+ array(pm8) = array(pm2)
229+ array(pm2) = temp
230+ endif
231+ if(compar(array(pm8), array(pm5)).eq.-1)then
232+ temp = array(pm8)
233+ array(pm8) = array(pm5)
234+ array(pm5) = temp
235+ endif
236+ endif
237+C Pivot assigned for medium and long length subarrays.
238+C Note that pm5 = mid
239+ piv = array(mid)
240+C Initialize pointers for partitioning
241+ i = lo-1
242+ j = hi+1
243+C Initialize counts of repeat values of pivot.
244+ a = 0
245+ d = 0
246+C Beginning of outer loop for placing pivot.
247+ 3 continue
248+C Scan up to find an element > piv.
249+ i = i + 1
250+C Check if pointers crossed.
251+ if(j.lt.i)goto 5
252+C Check if i pointer hit hi boundary.
253+ if(i.eq.hi)goto 4
254+C 
255+ if(compar(array(i), piv).eq.-1)goto 3
256+C Check for copies of pivot from scanning right. 
257+ if(compar(array(i), piv).eq.0)then
258+ array(i) = array(lo+a)
259+ array(lo+a) = piv
260+ a = a + 1
261+ goto 3 
262+ endif
263+C Beginning of innerloop for placing pivot.
264+ 4 continue
265+C Scan down to find an element < piv.
266+ j = j - 1
267+C Check if pointers crossed.
268+ if(j.lt.i)goto 5
269+ if(compar(array(j), piv).eq.1)goto 4
270+C Check for copies of pivot from scanning left. 
271+ if(compar(array(j), piv).eq.0)then
272+ array(j) = array(hi-d)
273+ array(hi-d) = piv
274+ d = d + 1
275+ goto 4 
276+ endif
277+C Check if pointers crossed.
278+ if(j.lt.i)goto 5
279+C Exchange elements
280+ temp = array(i)
281+ array(i) = array(j)
282+ array(j) = temp
283+C End of outermost loop for placing pivot.
284+ goto 3
285+C Insert all copies of pivot in appropriate place
286+ 5 s = MIN(a, j-lo-a+1)
287+ DO 6 k = 1, s
288+ array(lo-1+k) = array(i-k)
289+ array(i-k) = piv
290+ 6 CONTINUE
291+ s = MIN(d, hi-j-d)
292+ DO 7 k = 1, s
293+ array(hi+1-k) = array(j+k)
294+ array(j+k) = piv
295+ 7 CONTINUE
296+C Increase effective stack size
297+ tstack = tstack + 2 
298+C Push pointers to larger subarray on stack for later processing,
299+C process smaller subarray immediately. 
300+ if(tstack.gt.mstack) THEN
301+ WRITE(*,*)'Stack size is too small in quicksort fortran code quicksort_real_F77.F' 
302+ WRITE(*,*)'Are you sure you want to sort an array this long?'
303+ WRITE(*,*)'Your array has more than 2^63-1 entries?'
304+ WRITE(*,*)'If so, set mstack parameter to be at least:'
305+ WRITE(*,*)'2*ceiling(log_2(N))+2, for N = length of array,'
306+ WRITE(*,*)'and recompile this subroutine.'
307+ RETURN 
308+ endif
309+ if(hi-j-d-1.ge.j-a-lo)then
310+ bstack(tstack) = hi
311+ bstack(tstack-1) = MIN(j+d+1, hi)
312+ hi=MAX(j-a,lo)
313+ else
314+ bstack(tstack)=MAX(j-a,lo)
315+ bstack(tstack-1)=lo
316+ lo=MIN(j+d+1,hi)
317+ endif
318+C
319+C end of outermost if statement 
320+ endif
321+ goto 1
322+C END of subroutine quicksort
323+ END
tests/syntax-tests/source/Fortran (Fixed Form)/LICENSE.mdnew file+25
@@ -0,0 +1,25 @@
1+The `quicksort_real_F77.F` file has been added from https://github.com/jasonanema/Quicksort_Fortran77 under the following license:
2+
3+```text
4+MIT License
5+
6+Copyright (c) 2019 Jason Anema
7+
8+Permission is hereby granted, free of charge, to any person obtaining a copy
9+of this software and associated documentation files (the "Software"), to deal
10+in the Software without restriction, including without limitation the rights
11+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
12+copies of the Software, and to permit persons to whom the Software is
13+furnished to do so, subject to the following conditions:
14+
15+The above copyright notice and this permission notice shall be included in all
16+copies or substantial portions of the Software.
17+
18+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
19+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
20+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
21+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
22+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
23+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
24+SOFTWARE.
25+```
tests/syntax-tests/source/Fortran (Fixed Form)/quicksort_real_F77.Fnew file+323
@@ -0,0 +1,323 @@
1+C Fortran 77 implementation of a quicksort algorithm for arrays with
2+C real entries.
3+C ----------
4+C June 2019
5+C Jason Allen Anema, Ph.D.
6+C Division of Statistical Genomics
7+C Department of Genetics
8+C Washington University School of Medicine in St. Louis
9+C
10+C This work is partially supported NIH grant AG023746
11+C ----------
12+C Insertion sort is used for short arrays, as quicksort is slower on
13+C these.
14+C
15+C Hoare partition scheme is used (sweeping left and right), as it does
16+C three times fewer swaps on average that the Lamuto partition scheme.
17+C In conjunction with this, tripartite partition is performed
18+C concurrently (solving the "Dutch National Flag problem"). This avoids
19+C horrible runtimes on highly repetitive arrays. For example, without
20+C this, an array of random zeros and ones would have a runtime of
21+C O(N^2), but now has a runtime of O(N). The runtime for this algorthm
22+C on arrays with k highly repetitive entries is now O(kN).
23+C
24+C For medium length (sub)arrays, pivots are choosen using
25+C Median-of-Three, and those three items are sorted. For longer (sub)arrays
26+C the pseudomedian of nine (Median of medians). This avoids O(N^2) runtime on
27+C nonrandom inputs such as increasing and decreasing sequences.
28+C
29+C See Louis Bentley, Jon & McIlroy, Douglas. (1993). Engineering a Sort Function.
30+C Softw., Pract. Exper.. 23. 1249-1265. 10.1002/spe.4380231105 for details.
31+C
32+C The ordering on elements of the array are defined by a comparison
33+C function,compar, that is a user-supplied INTEGER*2 function of the form
34+C compar(a,b) which returns:
35+C -1 if a precedes b
36+C +1 if b precedes a
37+C 0 is a and b are considered equivalent
38+C and thus defines a total ordering.
39+C
40+C If one would like to use the standard order on integers, the
41+C compar function could be written in a file "compint.F" as:
42+C ----------------------------------------------------------------
43+C INTEGER*2 FUNCTION compint(a,b)
44+C INTEGER a, b
45+C if(a.lt.b)then
46+C compint = -1
47+C elseif(a.gt.b)then
48+C compint = +1
49+C else
50+C compint = 0
51+C endif
52+C END
53+C ----------------------------------------------------------------
54+C Then in your program, call quicksort with:
55+C call quicksort_real_F77(array, n, compint)
56+C
57+C The maximal length of an array in this implementation is (2^31-1),
58+C but can be changed to allow for length up to (2^63-1) by changing the
59+C data types of the relevant variables and constants. If you wish to
60+C sort longer arrays, of length N, you'll need to customize variable
61+C and constant types and set mstack to be at least (2*log_2(N)+2).
62+C
63+C ----------------------------------------------------------------
64+C Copyright 2019 Jason Allen Anema
65+C
66+C Permission is hereby granted, free of charge, to any person obtaining
67+C a copy of this software and associated documentation files (the "Software"),
68+C to deal in the Software without restriction, including without limitation the
69+C rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
70+C sell copies of the Software, and to permit persons to whom the Software is
71+C furnished to do so, subject to the following conditions:
72+C
73+C The above copyright notice and this permission notice shall be included
74+C in all copies or substantial portions of the Software.
75+C
76+C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
77+C OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
78+C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
79+C THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
80+C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
81+C FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
82+C IN THE SOFTWARE.
83+C -------------------------------------------------------------------
84+C
85+ SUBROUTINE quicksort_real_F77(array,n,compar)
86+ INTEGER n, maxins, maxmid, mstack
87+ REAL array(n)
88+ PARAMETER (maxins = 7, maxmid = 40, mstack = 128)
89+C maxins: maximal size of (sub)arrays to be sorted with
90+C insertion sort.
91+C maxmid: maximal size of (sub)arrays that will be quicksorted with
92+C Median-of-Three pivots.
93+C mstack: maximal size of required auxiliary storage (a stack), plus 2
94+C extra spots, which tracks the starts and ends of yet unsorted
95+C subarrays. mstack = 130 is large enough to handle arrays up to
96+C length 2^63-1. This maximal size follows from
97+C processing smaller arrays first and pigeonhole principal.
98+C
99+ INTEGER a, d, i, j, k, s, lo, mid, hi, tstack, bstack(mstack)
100+C a, d, i, j, k, s: indices
101+C lo, mid, and hi: their natural location in a (sub)array
102+C tstack: equal to twice the number of additional subarrays still
103+C needing to be sorted
104+C bstack: stack of the endpoints of unsorted subarrays
105+ INTEGER pm1, pm2, pm3, pm4, pm5, pm6, pm7, pm8, pm9
106+C for pseudomedian of nine positions in (sub)arrays
107+ REAL piv, temp
108+C piv is to store the pivot's value
109+C
110+ EXTERNAL compar
111+ INTEGER*2 compar
112+C compar is a user-supplied INTEGER*2 function of the form
113+C compar(a,b) which returns:
114+C -1 if a precedes b
115+C +1 if b precedes a
116+C 0 is a and b are considered equivalent
117+C and thus defines a total ordering.
118+ tstack = 0
119+ lo = 1
120+ hi = n
121+C
122+C Insertion sort subarrays of size maxins or less
123+ 1 if(hi-lo+1.le.maxins)then
124+ do 10, i = lo + 1, hi, 1
125+ temp = array(i)
126+ do 11 j = i - 1, lo, -1
127+ if(compar(array(j), temp).le.0)goto 2
128+ array(j+1)=array(j)
129+ 11 continue
130+ j = lo - 1
131+ 2 array(j+1) = temp
132+ 10 continue
133+ if(tstack.eq.0)return
134+C Pop the bstack, and start new partitioning
135+ hi = bstack(tstack)
136+ lo = bstack(tstack-1)
137+ tstack = tstack - 2
138+ else
139+C Use Median-of-Three as choice of pivot (median of lo, middle, hi)
140+C and reorder those elements appropriately when subarrays are medium
141+C length (between maxins and maxmid)
142+ mid = lo + (hi-lo)/2
143+ if(hi-lo.le.maxmid)then
144+ if(compar(array(mid), array(lo)).eq.-1)then
145+ temp = array(lo)
146+ array(lo) = array(mid)
147+ array(mid) = temp
148+ endif
149+ if(compar(array(hi), array(lo)).eq.-1)then
150+ temp = array(hi)
151+ array(hi) = array(lo)
152+ array(lo) = temp
153+ endif
154+ if(compar(array(hi), array(mid)).eq.-1)then
155+ temp = array(hi)
156+ array(hi) = array(mid)
157+ array(mid) = temp
158+ endif
159+C Use pseudomedian of nine (Median of medians) as choice of pivot when
160+C subarrays are longer than maxmid. Note that doing it this way requires only 12
161+C comparisons for finding the pivot.
162+ elseif(hi-lo+1.gt.maxmid)then
163+ pm1 = lo
164+ pm5 = lo + (hi-lo)/2
165+ pm9 = hi
166+ pm3 = lo + (pm5-lo)/2
167+ pm7 = pm5 + (hi-pm5)/2
168+ pm2 = lo + (pm3-lo)/2
169+ pm4 = pm3 + (pm5-pm3)/2
170+ pm6 = pm5 + (pm7-pm5)/2
171+ pm8 = pm7 + (pm9-pm7)/2
172+C Median and sorting for pm1, pm2, pm3
173+ if(compar(array(pm2), array(pm1)).eq.-1)then
174+ temp = array(pm1)
175+ array(pm1) = array(pm2)
176+ array(pm2) = temp
177+ endif
178+ if(compar(array(pm3), array(pm1)).eq.-1)then
179+ temp = array(pm3)
180+ array(pm3) = array(pm1)
181+ array(pm1) = temp
182+ endif
183+ if(compar(array(pm3), array(pm2)).eq.-1)then
184+ temp = array(pm3)
185+ array(pm3) = array(pm2)
186+ array(pm2) = temp
187+ endif
188+C Median and sorting for pm4, pm5, pm6
189+ if(compar(array(pm5), array(pm4)).eq.-1)then
190+ temp = array(pm4)
191+ array(pm4) = array(pm5)
192+ array(pm5) = temp
193+ endif
194+ if(compar(array(pm6), array(pm4)).eq.-1)then
195+ temp = array(pm6)
196+ array(pm6) = array(pm4)
197+ array(pm4) = temp
198+ endif
199+ if(compar(array(pm6), array(pm5)).eq.-1)then
200+ temp = array(pm6)
201+ array(pm6) = array(pm5)
202+ array(pm5) = temp
203+ endif
204+C Median and sorting for pm7, pm8, pm9
205+ if(compar(array(pm8), array(pm7)).eq.-1)then
206+ temp = array(pm7)
207+ array(pm7) = array(pm8)
208+ array(pm8) = temp
209+ endif
210+ if(compar(array(pm9), array(pm7)).eq.-1)then
211+ temp = array(pm9)
212+ array(pm9) = array(pm7)
213+ array(pm7) = temp
214+ endif
215+ if(compar(array(pm9), array(pm8)).eq.-1)then
216+ temp = array(pm9)
217+ array(pm9) = array(pm8)
218+ array(pm8) = temp
219+ endif
220+C Median of the medians (which are now pm2, pm5, pm8)
221+ if(compar(array(pm5), array(pm2)).eq.-1)then
222+ temp = array(pm2)
223+ array(pm2) = array(pm5)
224+ array(pm5) = temp
225+ endif
226+ if(compar(array(pm8), array(pm2)).eq.-1)then
227+ temp = array(pm8)
228+ array(pm8) = array(pm2)
229+ array(pm2) = temp
230+ endif
231+ if(compar(array(pm8), array(pm5)).eq.-1)then
232+ temp = array(pm8)
233+ array(pm8) = array(pm5)
234+ array(pm5) = temp
235+ endif
236+ endif
237+C Pivot assigned for medium and long length subarrays.
238+C Note that pm5 = mid
239+ piv = array(mid)
240+C Initialize pointers for partitioning
241+ i = lo-1
242+ j = hi+1
243+C Initialize counts of repeat values of pivot.
244+ a = 0
245+ d = 0
246+C Beginning of outer loop for placing pivot.
247+ 3 continue
248+C Scan up to find an element > piv.
249+ i = i + 1
250+C Check if pointers crossed.
251+ if(j.lt.i)goto 5
252+C Check if i pointer hit hi boundary.
253+ if(i.eq.hi)goto 4
254+C
255+ if(compar(array(i), piv).eq.-1)goto 3
256+C Check for copies of pivot from scanning right.
257+ if(compar(array(i), piv).eq.0)then
258+ array(i) = array(lo+a)
259+ array(lo+a) = piv
260+ a = a + 1
261+ goto 3
262+ endif
263+C Beginning of innerloop for placing pivot.
264+ 4 continue
265+C Scan down to find an element < piv.
266+ j = j - 1
267+C Check if pointers crossed.
268+ if(j.lt.i)goto 5
269+ if(compar(array(j), piv).eq.1)goto 4
270+C Check for copies of pivot from scanning left.
271+ if(compar(array(j), piv).eq.0)then
272+ array(j) = array(hi-d)
273+ array(hi-d) = piv
274+ d = d + 1
275+ goto 4
276+ endif
277+C Check if pointers crossed.
278+ if(j.lt.i)goto 5
279+C Exchange elements
280+ temp = array(i)
281+ array(i) = array(j)
282+ array(j) = temp
283+C End of outermost loop for placing pivot.
284+ goto 3
285+C Insert all copies of pivot in appropriate place
286+ 5 s = MIN(a, j-lo-a+1)
287+ DO 6 k = 1, s
288+ array(lo-1+k) = array(i-k)
289+ array(i-k) = piv
290+ 6 CONTINUE
291+ s = MIN(d, hi-j-d)
292+ DO 7 k = 1, s
293+ array(hi+1-k) = array(j+k)
294+ array(j+k) = piv
295+ 7 CONTINUE
296+C Increase effective stack size
297+ tstack = tstack + 2
298+C Push pointers to larger subarray on stack for later processing,
299+C process smaller subarray immediately.
300+ if(tstack.gt.mstack) THEN
301+ WRITE(*,*)'Stack size is too small in quicksort fortran code quicksort_real_F77.F'
302+ WRITE(*,*)'Are you sure you want to sort an array this long?'
303+ WRITE(*,*)'Your array has more than 2^63-1 entries?'
304+ WRITE(*,*)'If so, set mstack parameter to be at least:'
305+ WRITE(*,*)'2*ceiling(log_2(N))+2, for N = length of array,'
306+ WRITE(*,*)'and recompile this subroutine.'
307+ RETURN
308+ endif
309+ if(hi-j-d-1.ge.j-a-lo)then
310+ bstack(tstack) = hi
311+ bstack(tstack-1) = MIN(j+d+1, hi)
312+ hi=MAX(j-a,lo)
313+ else
314+ bstack(tstack)=MAX(j-a,lo)
315+ bstack(tstack-1)=lo
316+ lo=MIN(j+d+1,hi)
317+ endif
318+C
319+C end of outermost if statement
320+ endif
321+ goto 1
322+C END of subroutine quicksort
323+ END