Add Fortran (Fixed Form) syntax test file

tests/syntax-tests/highlighted/Fortran (Fixed Form)/quicksort_real_F77.F

@@ -0,0 +1,323 @@
+C  Fortran 77 implementation of a quicksort algorithm for arrays with
+C  real entries.
+C  ----------
+C  June 2019 
+C  Jason Allen Anema, Ph.D.
+C  Division of Statistical Genomics
+C  Department of Genetics
+C  Washington University School of Medicine in St. Louis
+C 
+C  This work is partially supported NIH grant AG023746
+C  ----------
+C  Insertion sort is used for short arrays, as quicksort is slower on
+C  these.
+C  
+C  Hoare partition scheme is used (sweeping left and right), as it does
+C  three times fewer swaps on average that the Lamuto partition scheme.
+C  In conjunction with this, tripartite partition is performed
+C  concurrently (solving the "Dutch National Flag problem"). This avoids 
+C  horrible runtimes on highly repetitive arrays. For example, without  
+C  this, an array of random zeros and ones would have a runtime of
+C  O(N^2), but now has a runtime of O(N). The runtime for this algorthm
+C  on arrays with k highly repetitive entries is now O(kN).
+C    
+C  For medium length (sub)arrays, pivots are choosen using
+C  Median-of-Three, and those three items are sorted. For longer (sub)arrays
+C  the pseudomedian of nine (Median of medians). This avoids O(N^2) runtime on
+C  nonrandom inputs such as increasing and decreasing sequences. 
+C
+C  See Louis Bentley, Jon & McIlroy, Douglas. (1993). Engineering a Sort Function.
+C  Softw., Pract. Exper.. 23. 1249-1265. 10.1002/spe.4380231105 for details. 
+C
+C  The ordering on elements of the array are defined by a comparison
+C  function,compar, that is a user-supplied INTEGER*2 function of the form
+C           compar(a,b) which returns:
+C              -1 if a precedes b
+C              +1 if b precedes a
+C              0 is a and b are considered equivalent
+C  and thus defines a total ordering.
+C  
+C  If one would like to use the standard order on integers, the
+C  compar function could be written in a file "compint.F" as:
+C  ----------------------------------------------------------------
+C      INTEGER*2 FUNCTION compint(a,b)
+C      INTEGER a, b
+C      if(a.lt.b)then
+C         compint = -1
+C      elseif(a.gt.b)then
+C         compint = +1
+C      else
+C         compint = 0
+C      endif
+C      END
+C  ----------------------------------------------------------------
+C  Then in your program, call quicksort with:
+C      call quicksort_real_F77(array, n, compint)
+C
+C  The maximal length of an array in this implementation is (2^31-1),
+C  but can be changed to allow for length up to (2^63-1) by changing the
+C  data types of the relevant variables and constants. If you wish to 
+C  sort longer arrays, of length N, you'll need to customize variable 
+C  and constant types and set mstack to be at least (2*log_2(N)+2). 
+C
+C  ----------------------------------------------------------------
+C  Copyright 2019 Jason Allen Anema
+C  
+C  Permission is hereby granted, free of charge, to any person obtaining
+C  a copy of this software and associated documentation files (the "Software"),
+C  to deal in the Software without restriction, including without limitation the
+C  rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+C  sell copies of the Software, and to permit persons to whom the Software is
+C  furnished to do so, subject to the following conditions:
+C
+C  The above copyright notice and this permission notice shall be included
+C  in all copies or substantial portions of the Software.
+C
+C  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+C  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+C  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+C  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+C  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+C  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+C  IN THE SOFTWARE.
+C  -------------------------------------------------------------------
+C
+      SUBROUTINE quicksort_real_F77(array,n,compar)
+      INTEGER n, maxins, maxmid, mstack
+      REAL array(n)
+      PARAMETER (maxins = 7, maxmid = 40, mstack = 128)
+C        maxins: maximal size of (sub)arrays to be sorted with
+C           insertion sort.
+C        maxmid: maximal size of (sub)arrays that will be quicksorted with
+C           Median-of-Three pivots.
+C        mstack: maximal size of required auxiliary storage (a stack), plus 2 
+C           extra spots, which tracks the starts and ends of yet unsorted 
+C           subarrays. mstack = 130 is large enough to handle arrays up to 
+C           length 2^63-1. This maximal size follows from
+C           processing smaller arrays first and pigeonhole principal.
+C 
+      INTEGER  a, d, i, j, k, s, lo, mid, hi, tstack, bstack(mstack)
+C        a, d, i, j, k, s: indices
+C        lo, mid, and hi: their natural location in a (sub)array
+C        tstack:  equal to twice the number of additional subarrays still 
+C          needing to be sorted
+C        bstack: stack of the endpoints of unsorted subarrays
+      INTEGER pm1, pm2, pm3, pm4, pm5, pm6, pm7, pm8, pm9
+C        for pseudomedian of nine positions in (sub)arrays
+      REAL piv, temp
+C        piv is to store the pivot's value
+C    
+      EXTERNAL compar
+      INTEGER*2 compar
+C        compar is a user-supplied INTEGER*2 function of the form
+C           compar(a,b) which returns:
+C              -1 if a precedes b
+C              +1 if b precedes a
+C              0 is a and b are considered equivalent
+C           and thus defines a total ordering. 
+      tstack = 0
+      lo = 1 
+      hi = n
+C
+C  Insertion sort subarrays of size maxins or less
+ 1    if(hi-lo+1.le.maxins)then
+         do 10, i = lo + 1, hi, 1
+            temp = array(i)
+            do 11 j = i - 1, lo, -1
+               if(compar(array(j), temp).le.0)goto 2
+               array(j+1)=array(j)
+ 11         continue
+            j = lo - 1
+ 2          array(j+1) = temp
+ 10      continue
+         if(tstack.eq.0)return
+C  Pop the bstack, and start new partitioning
+         hi = bstack(tstack)
+         lo = bstack(tstack-1)
+         tstack = tstack - 2
+      else
+C  Use Median-of-Three as choice of pivot (median of lo, middle, hi)
+C  and reorder those elements appropriately when subarrays are medium
+C  length (between maxins and maxmid)
+         mid = lo +  (hi-lo)/2
+         if(hi-lo.le.maxmid)then
+            if(compar(array(mid), array(lo)).eq.-1)then
+               temp = array(lo)
+               array(lo) = array(mid)
+               array(mid) = temp
+            endif
+            if(compar(array(hi), array(lo)).eq.-1)then
+               temp = array(hi)
+               array(hi) = array(lo)
+               array(lo) = temp
+            endif
+            if(compar(array(hi), array(mid)).eq.-1)then
+               temp = array(hi)
+               array(hi) = array(mid)
+               array(mid) = temp
+            endif
+C  Use pseudomedian of nine (Median of medians) as choice of pivot when 
+C  subarrays are longer than maxmid. Note that doing it this way requires only 12
+C  comparisons for finding the pivot.
+         elseif(hi-lo+1.gt.maxmid)then
+            pm1 = lo
+            pm5 = lo + (hi-lo)/2
+            pm9 = hi
+            pm3 = lo + (pm5-lo)/2
+            pm7 = pm5 + (hi-pm5)/2
+            pm2 = lo + (pm3-lo)/2
+            pm4 = pm3 + (pm5-pm3)/2
+            pm6 = pm5 + (pm7-pm5)/2
+            pm8 = pm7 + (pm9-pm7)/2
+C  Median and sorting for pm1, pm2, pm3
+            if(compar(array(pm2), array(pm1)).eq.-1)then
+               temp = array(pm1)
+               array(pm1) = array(pm2)
+               array(pm2) = temp
+            endif
+            if(compar(array(pm3), array(pm1)).eq.-1)then
+               temp = array(pm3)
+               array(pm3) = array(pm1)
+               array(pm1) = temp
+            endif
+            if(compar(array(pm3), array(pm2)).eq.-1)then
+               temp = array(pm3)
+               array(pm3) = array(pm2)
+               array(pm2) = temp
+            endif
+C  Median and sorting for pm4, pm5, pm6
+            if(compar(array(pm5), array(pm4)).eq.-1)then
+               temp = array(pm4)
+               array(pm4) = array(pm5)
+               array(pm5) = temp
+            endif
+            if(compar(array(pm6), array(pm4)).eq.-1)then
+               temp = array(pm6)
+               array(pm6) = array(pm4)
+               array(pm4) = temp
+            endif
+            if(compar(array(pm6), array(pm5)).eq.-1)then
+               temp = array(pm6)
+               array(pm6) = array(pm5)
+               array(pm5) = temp
+            endif
+C  Median and sorting for pm7, pm8, pm9
+            if(compar(array(pm8), array(pm7)).eq.-1)then
+               temp = array(pm7)
+               array(pm7) = array(pm8)
+               array(pm8) = temp
+            endif
+            if(compar(array(pm9), array(pm7)).eq.-1)then
+               temp = array(pm9)
+               array(pm9) = array(pm7)
+               array(pm7) = temp
+            endif
+            if(compar(array(pm9), array(pm8)).eq.-1)then
+               temp = array(pm9)
+               array(pm9) = array(pm8)
+               array(pm8) = temp
+            endif
+C Median of the medians (which are now pm2, pm5, pm8)
+            if(compar(array(pm5), array(pm2)).eq.-1)then
+               temp = array(pm2)
+               array(pm2) = array(pm5)
+               array(pm5) = temp
+            endif
+            if(compar(array(pm8), array(pm2)).eq.-1)then
+               temp = array(pm8)
+               array(pm8) = array(pm2)
+               array(pm2) = temp
+            endif
+            if(compar(array(pm8), array(pm5)).eq.-1)then
+               temp = array(pm8)
+               array(pm8) = array(pm5)
+               array(pm5) = temp
+            endif
+         endif
+C  Pivot assigned for medium and long length subarrays.
+C  Note that pm5 = mid
+             piv = array(mid)
+C  Initialize pointers for partitioning
+            i = lo-1
+            j = hi+1
+C  Initialize counts of repeat values of pivot.
+            a = 0
+            d = 0
+C  Beginning of outer loop for placing pivot.
+ 3       continue
+C  Scan up to find an element > piv.
+            i = i + 1
+C  Check if pointers crossed.
+         if(j.lt.i)goto 5
+C  Check if i pointer hit hi boundary.
+         if(i.eq.hi)goto 4
+C      
+         if(compar(array(i), piv).eq.-1)goto 3
+C Check for copies of pivot from scanning right. 
+         if(compar(array(i), piv).eq.0)then
+             array(i) = array(lo+a)
+             array(lo+a) = piv
+             a = a + 1
+             goto 3 
+         endif
+C  Beginning of innerloop for placing pivot.
+ 4       continue
+C  Scan down to find an element < piv.
+            j = j - 1
+C  Check if pointers crossed.
+         if(j.lt.i)goto 5
+         if(compar(array(j), piv).eq.1)goto 4
+C  Check for copies of pivot from scanning left. 
+         if(compar(array(j), piv).eq.0)then
+            array(j) = array(hi-d)
+            array(hi-d) = piv
+            d = d + 1
+            goto 4      
+         endif
+C Check if pointers crossed.
+         if(j.lt.i)goto 5
+C  Exchange elements
+         temp = array(i)
+         array(i) = array(j)
+         array(j) = temp
+C  End of outermost loop for placing pivot.
+         goto 3
+C  Insert all copies of pivot in appropriate place
+ 5       s = MIN(a, j-lo-a+1)
+         DO 6 k = 1, s
+            array(lo-1+k) = array(i-k)
+            array(i-k) = piv
+ 6       CONTINUE
+         s = MIN(d, hi-j-d)
+         DO 7 k = 1, s
+            array(hi+1-k) = array(j+k)
+            array(j+k) = piv
+ 7       CONTINUE
+C  Increase effective stack size
+         tstack = tstack + 2 
+C  Push pointers to larger subarray on stack for later processing,
+C  process smaller subarray immediately. 
+         if(tstack.gt.mstack) THEN
+           WRITE(*,*)'Stack size is too small in quicksort fortran code quicksort_real_F77.F' 
+           WRITE(*,*)'Are you sure you want to sort an array this long?'
+           WRITE(*,*)'Your array has more than 2^63-1  entries?'
+           WRITE(*,*)'If so, set mstack parameter to be at least:'
+           WRITE(*,*)'2*ceiling(log_2(N))+2, for N = length of array,'
+           WRITE(*,*)'and recompile this subroutine.'
+           RETURN 
+         endif
+         if(hi-j-d-1.ge.j-a-lo)then
+            bstack(tstack) = hi
+            bstack(tstack-1) = MIN(j+d+1, hi)
+            hi=MAX(j-a,lo)
+         else
+            bstack(tstack)=MAX(j-a,lo)
+            bstack(tstack-1)=lo
+            lo=MIN(j+d+1,hi)
+         endif
+C
+C     end of outermost if statement 
+      endif
+      goto 1
+C     END of subroutine quicksort
+      END

tests/syntax-tests/source/Fortran (Fixed Form)/LICENSE.md

@@ -0,0 +1,25 @@
+The `quicksort_real_F77.F` file has been added from https://github.com/jasonanema/Quicksort_Fortran77 under the following license:
+
+```text
+MIT License
+
+Copyright (c) 2019 Jason Anema
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+```

tests/syntax-tests/source/Fortran (Fixed Form)/quicksort_real_F77.F

@@ -0,0 +1,323 @@
+C  Fortran 77 implementation of a quicksort algorithm for arrays with
+C  real entries.
+C  ----------
+C  June 2019 
+C  Jason Allen Anema, Ph.D.
+C  Division of Statistical Genomics
+C  Department of Genetics
+C  Washington University School of Medicine in St. Louis
+C 
+C  This work is partially supported NIH grant AG023746
+C  ----------
+C  Insertion sort is used for short arrays, as quicksort is slower on
+C  these.
+C  
+C  Hoare partition scheme is used (sweeping left and right), as it does
+C  three times fewer swaps on average that the Lamuto partition scheme.
+C  In conjunction with this, tripartite partition is performed
+C  concurrently (solving the "Dutch National Flag problem"). This avoids 
+C  horrible runtimes on highly repetitive arrays. For example, without  
+C  this, an array of random zeros and ones would have a runtime of
+C  O(N^2), but now has a runtime of O(N). The runtime for this algorthm
+C  on arrays with k highly repetitive entries is now O(kN).
+C    
+C  For medium length (sub)arrays, pivots are choosen using
+C  Median-of-Three, and those three items are sorted. For longer (sub)arrays
+C  the pseudomedian of nine (Median of medians). This avoids O(N^2) runtime on
+C  nonrandom inputs such as increasing and decreasing sequences. 
+C
+C  See Louis Bentley, Jon & McIlroy, Douglas. (1993). Engineering a Sort Function.
+C  Softw., Pract. Exper.. 23. 1249-1265. 10.1002/spe.4380231105 for details. 
+C
+C  The ordering on elements of the array are defined by a comparison
+C  function,compar, that is a user-supplied INTEGER*2 function of the form
+C           compar(a,b) which returns:
+C              -1 if a precedes b
+C              +1 if b precedes a
+C              0 is a and b are considered equivalent
+C  and thus defines a total ordering.
+C  
+C  If one would like to use the standard order on integers, the
+C  compar function could be written in a file "compint.F" as:
+C  ----------------------------------------------------------------
+C      INTEGER*2 FUNCTION compint(a,b)
+C      INTEGER a, b
+C      if(a.lt.b)then
+C         compint = -1
+C      elseif(a.gt.b)then
+C         compint = +1
+C      else
+C         compint = 0
+C      endif
+C      END
+C  ----------------------------------------------------------------
+C  Then in your program, call quicksort with:
+C      call quicksort_real_F77(array, n, compint)
+C
+C  The maximal length of an array in this implementation is (2^31-1),
+C  but can be changed to allow for length up to (2^63-1) by changing the
+C  data types of the relevant variables and constants. If you wish to 
+C  sort longer arrays, of length N, you'll need to customize variable 
+C  and constant types and set mstack to be at least (2*log_2(N)+2). 
+C
+C  ----------------------------------------------------------------
+C  Copyright 2019 Jason Allen Anema
+C  
+C  Permission is hereby granted, free of charge, to any person obtaining
+C  a copy of this software and associated documentation files (the "Software"),
+C  to deal in the Software without restriction, including without limitation the
+C  rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
+C  sell copies of the Software, and to permit persons to whom the Software is
+C  furnished to do so, subject to the following conditions:
+C
+C  The above copyright notice and this permission notice shall be included
+C  in all copies or substantial portions of the Software.
+C
+C  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+C  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+C  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+C  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+C  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+C  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
+C  IN THE SOFTWARE.
+C  -------------------------------------------------------------------
+C
+      SUBROUTINE quicksort_real_F77(array,n,compar)
+      INTEGER n, maxins, maxmid, mstack
+      REAL array(n)
+      PARAMETER (maxins = 7, maxmid = 40, mstack = 128)
+C        maxins: maximal size of (sub)arrays to be sorted with
+C           insertion sort.
+C        maxmid: maximal size of (sub)arrays that will be quicksorted with
+C           Median-of-Three pivots.
+C        mstack: maximal size of required auxiliary storage (a stack), plus 2 
+C           extra spots, which tracks the starts and ends of yet unsorted 
+C           subarrays. mstack = 130 is large enough to handle arrays up to 
+C           length 2^63-1. This maximal size follows from
+C           processing smaller arrays first and pigeonhole principal.
+C 
+      INTEGER  a, d, i, j, k, s, lo, mid, hi, tstack, bstack(mstack)
+C        a, d, i, j, k, s: indices
+C        lo, mid, and hi: their natural location in a (sub)array
+C        tstack:  equal to twice the number of additional subarrays still 
+C          needing to be sorted
+C        bstack: stack of the endpoints of unsorted subarrays
+      INTEGER pm1, pm2, pm3, pm4, pm5, pm6, pm7, pm8, pm9
+C        for pseudomedian of nine positions in (sub)arrays
+      REAL piv, temp
+C        piv is to store the pivot's value
+C    
+      EXTERNAL compar
+      INTEGER*2 compar
+C        compar is a user-supplied INTEGER*2 function of the form
+C           compar(a,b) which returns:
+C              -1 if a precedes b
+C              +1 if b precedes a
+C              0 is a and b are considered equivalent
+C           and thus defines a total ordering. 
+      tstack = 0
+      lo = 1 
+      hi = n
+C
+C  Insertion sort subarrays of size maxins or less
+ 1    if(hi-lo+1.le.maxins)then
+         do 10, i = lo + 1, hi, 1
+            temp = array(i)
+            do 11 j = i - 1, lo, -1
+               if(compar(array(j), temp).le.0)goto 2
+               array(j+1)=array(j)
+ 11         continue
+            j = lo - 1
+ 2          array(j+1) = temp
+ 10      continue
+         if(tstack.eq.0)return
+C  Pop the bstack, and start new partitioning
+         hi = bstack(tstack)
+         lo = bstack(tstack-1)
+         tstack = tstack - 2
+      else
+C  Use Median-of-Three as choice of pivot (median of lo, middle, hi)
+C  and reorder those elements appropriately when subarrays are medium
+C  length (between maxins and maxmid)
+         mid = lo +  (hi-lo)/2
+         if(hi-lo.le.maxmid)then
+            if(compar(array(mid), array(lo)).eq.-1)then
+               temp = array(lo)
+               array(lo) = array(mid)
+               array(mid) = temp
+            endif
+            if(compar(array(hi), array(lo)).eq.-1)then
+               temp = array(hi)
+               array(hi) = array(lo)
+               array(lo) = temp
+            endif
+            if(compar(array(hi), array(mid)).eq.-1)then
+               temp = array(hi)
+               array(hi) = array(mid)
+               array(mid) = temp
+            endif
+C  Use pseudomedian of nine (Median of medians) as choice of pivot when 
+C  subarrays are longer than maxmid. Note that doing it this way requires only 12
+C  comparisons for finding the pivot.
+         elseif(hi-lo+1.gt.maxmid)then
+            pm1 = lo
+            pm5 = lo + (hi-lo)/2
+            pm9 = hi
+            pm3 = lo + (pm5-lo)/2
+            pm7 = pm5 + (hi-pm5)/2
+            pm2 = lo + (pm3-lo)/2
+            pm4 = pm3 + (pm5-pm3)/2
+            pm6 = pm5 + (pm7-pm5)/2
+            pm8 = pm7 + (pm9-pm7)/2
+C  Median and sorting for pm1, pm2, pm3
+            if(compar(array(pm2), array(pm1)).eq.-1)then
+               temp = array(pm1)
+               array(pm1) = array(pm2)
+               array(pm2) = temp
+            endif
+            if(compar(array(pm3), array(pm1)).eq.-1)then
+               temp = array(pm3)
+               array(pm3) = array(pm1)
+               array(pm1) = temp
+            endif
+            if(compar(array(pm3), array(pm2)).eq.-1)then
+               temp = array(pm3)
+               array(pm3) = array(pm2)
+               array(pm2) = temp
+            endif
+C  Median and sorting for pm4, pm5, pm6
+            if(compar(array(pm5), array(pm4)).eq.-1)then
+               temp = array(pm4)
+               array(pm4) = array(pm5)
+               array(pm5) = temp
+            endif
+            if(compar(array(pm6), array(pm4)).eq.-1)then
+               temp = array(pm6)
+               array(pm6) = array(pm4)
+               array(pm4) = temp
+            endif
+            if(compar(array(pm6), array(pm5)).eq.-1)then
+               temp = array(pm6)
+               array(pm6) = array(pm5)
+               array(pm5) = temp
+            endif
+C  Median and sorting for pm7, pm8, pm9
+            if(compar(array(pm8), array(pm7)).eq.-1)then
+               temp = array(pm7)
+               array(pm7) = array(pm8)
+               array(pm8) = temp
+            endif
+            if(compar(array(pm9), array(pm7)).eq.-1)then
+               temp = array(pm9)
+               array(pm9) = array(pm7)
+               array(pm7) = temp
+            endif
+            if(compar(array(pm9), array(pm8)).eq.-1)then
+               temp = array(pm9)
+               array(pm9) = array(pm8)
+               array(pm8) = temp
+            endif
+C Median of the medians (which are now pm2, pm5, pm8)
+            if(compar(array(pm5), array(pm2)).eq.-1)then
+               temp = array(pm2)
+               array(pm2) = array(pm5)
+               array(pm5) = temp
+            endif
+            if(compar(array(pm8), array(pm2)).eq.-1)then
+               temp = array(pm8)
+               array(pm8) = array(pm2)
+               array(pm2) = temp
+            endif
+            if(compar(array(pm8), array(pm5)).eq.-1)then
+               temp = array(pm8)
+               array(pm8) = array(pm5)
+               array(pm5) = temp
+            endif
+         endif
+C  Pivot assigned for medium and long length subarrays.
+C  Note that pm5 = mid
+             piv = array(mid)
+C  Initialize pointers for partitioning
+            i = lo-1
+            j = hi+1
+C  Initialize counts of repeat values of pivot.
+            a = 0
+            d = 0
+C  Beginning of outer loop for placing pivot.
+ 3       continue
+C  Scan up to find an element > piv.
+            i = i + 1
+C  Check if pointers crossed.
+         if(j.lt.i)goto 5
+C  Check if i pointer hit hi boundary.
+         if(i.eq.hi)goto 4
+C      
+         if(compar(array(i), piv).eq.-1)goto 3
+C Check for copies of pivot from scanning right. 
+         if(compar(array(i), piv).eq.0)then
+             array(i) = array(lo+a)
+             array(lo+a) = piv
+             a = a + 1
+             goto 3 
+         endif
+C  Beginning of innerloop for placing pivot.
+ 4       continue
+C  Scan down to find an element < piv.
+            j = j - 1
+C  Check if pointers crossed.
+         if(j.lt.i)goto 5
+         if(compar(array(j), piv).eq.1)goto 4
+C  Check for copies of pivot from scanning left. 
+         if(compar(array(j), piv).eq.0)then
+            array(j) = array(hi-d)
+            array(hi-d) = piv
+            d = d + 1
+            goto 4      
+         endif
+C Check if pointers crossed.
+         if(j.lt.i)goto 5
+C  Exchange elements
+         temp = array(i)
+         array(i) = array(j)
+         array(j) = temp
+C  End of outermost loop for placing pivot.
+         goto 3
+C  Insert all copies of pivot in appropriate place
+ 5       s = MIN(a, j-lo-a+1)
+         DO 6 k = 1, s
+            array(lo-1+k) = array(i-k)
+            array(i-k) = piv
+ 6       CONTINUE
+         s = MIN(d, hi-j-d)
+         DO 7 k = 1, s
+            array(hi+1-k) = array(j+k)
+            array(j+k) = piv
+ 7       CONTINUE
+C  Increase effective stack size
+         tstack = tstack + 2 
+C  Push pointers to larger subarray on stack for later processing,
+C  process smaller subarray immediately. 
+         if(tstack.gt.mstack) THEN
+           WRITE(*,*)'Stack size is too small in quicksort fortran code quicksort_real_F77.F' 
+           WRITE(*,*)'Are you sure you want to sort an array this long?'
+           WRITE(*,*)'Your array has more than 2^63-1  entries?'
+           WRITE(*,*)'If so, set mstack parameter to be at least:'
+           WRITE(*,*)'2*ceiling(log_2(N))+2, for N = length of array,'
+           WRITE(*,*)'and recompile this subroutine.'
+           RETURN 
+         endif
+         if(hi-j-d-1.ge.j-a-lo)then
+            bstack(tstack) = hi
+            bstack(tstack-1) = MIN(j+d+1, hi)
+            hi=MAX(j-a,lo)
+         else
+            bstack(tstack)=MAX(j-a,lo)
+            bstack(tstack-1)=lo
+            lo=MIN(j+d+1,hi)
+         endif
+C
+C     end of outermost if statement 
+      endif
+      goto 1
+C     END of subroutine quicksort
+      END