c	From nwang@PICARD.tamu.eduThu Jul 27 17:27:22 1995
c	Date: Thu, 27 Jul 95 17:28:44 CDT
c	From: Naisyin Wang <nwang@PICARD.tamu.edu>
c	To: sengupta@eesun1.tamu.edu
c	Subject: qsort.f

c...	Modified by Manidip Sengupta, as follows:
c...	[1] No implicit variables
c...	[2] Sorts a double prec. array

c-------- --------- --------- --------- --------- --------- --------- --|
       subroutine dqksort(x,n)
c*********************************************************************
c
c   Quicksort a real array x of length n.
c
c********************************************************************
c
c      dimension x(n),ll(20),lr(20)

	implicit none
	integer i, j, ns, nl, nm, n, nleft, nright, itemp, ll(20),lr(20)
	double precision ax, x(n)

c
c   ns is # of pieces left to partition
c
c   ll and lr are left and right borders of pieces
c
       ns=1
       ll(1)=1
       lr(1)=n
c
c   After splitting we come here. Stop when no pieces left.
c
   5   if(ns.eq.0) go to 99
c
       i=ll(ns)
       j=lr(ns)
       nl=j-i+1
c
c   Sort pieces of size 1 or 2:
c
       if(nl.le.2) then
              ns=ns-1
              if(nl.eq.2.and.x(i).gt.x(j)) call dswap(x(i),x(j))
              go to 5
       endif
c
c   Sort first, ``middle'' and last elements:
c
       nm=(i+j)/2
       if(x(i).gt.x(nm)) call dswap(x(i),x(nm))
       if(x(nm).gt.x(j)) call dswap(x(nm),x(j))
       if(x(i).gt.x(nm)) call dswap(x(i),x(nm))
c
c   If piece of size 3, it's now sorted
c
       if(nl.eq.3) then
              ns=ns-1
              go to 5
       endif
c
c   Put middle (target) into 1st element and keep a copy for comparisons
c
       ax=x(nm)
       call dswap(x(i),x(nm))
c
c   Now we'll look for 1st one from left (starting with 2nd)  > target 
c   and first from right < target. If no such pair we end up at 20.
c
       itemp=i
       i=i+1
c
  10   if(x(i).le.ax) then
              if(i.eq.j) go to 20
              i=i+1
              go to 10
       endif
c
  15   if(x(j).ge.ax) then
              if(i.eq.j) go to 20
              j=j-1
              go to 15
       endif
c
c   Can only get to here if we found a pair to swap
c
       call dswap(x(i),x(j))
c
c   Only go back to look for another pair if there's a chance to get one
c
       if(j-i.le.1) then
              go to 20
       else
              i=i+1
              j=j-1
              go to 10
       endif
c
c   Get to here when no more pairs to swap. i might be 1 space too far
c   to right 
c
  20   if(x(i).gt.ax) i=i-1
       call dswap(x(itemp),x(i))
c
c   Now put shorter piece at end of ll and lr (this guarantees small
c   number of pieces at any one time). A piece might be empty, but that
c   is taken care of above.
c
       nleft=i-ll(ns)
       nright=lr(ns)-i
       ns=ns+1
       if(nleft.le.nright) then
              ll(ns)=ll(ns-1)
              ll(ns-1)=i+1
              lr(ns)=i-1
       else
              lr(ns)=lr(ns-1)
              lr(ns-1)=i-1
              ll(ns)=i+1
       endif
       go to 5
c
  99   continue
       return
       end


       subroutine dswap(a,b)
c*********************************************************************
c
c   Swap two double precision numbers.
c
c********************************************************************
c
	double precision a, b, c
       c=b
       b=a
       a=c
c
       return
       end