subroutine quick_sort(n1, n2, ld, exch, comp, userarg)

C Quicksort Algorithm        M. Zolliker Dec. 89
C
C N = number of elements
C LD, EXCH and COMP are subroutines
C
C SUBROUTINE LD(I,userarg):   load element no. I into a buffer
C
C SUBROUTINE EXCH(I,J,userarg): exchange elements no. I and J
C
C SUBROUTINE COMP(I, TST,userarg): compare element no. I with buffer, result is TST
C
C	    <                 <
C	TST = 0  if element I = buffer
C	    >                 >
C
C	for example for comparing reals:
C
C	subroutine comp(i, tst)
C	common buffer
C	if (element(i) .ne. buffer) then
C	  if (element(i) .gt. buffer) then
C           tst=1
C	  else
C           tst=-1
C	  endif
C	else
C	  tst=0
C	endif
C	end
C
	implicit none
	integer nstack
	parameter (nstack=30)	! sequence stack dim. for 2**30 elements
	integer n1, n2, i, j, k, l, m, tst
	integer kst(nstack), lst(nstack), nst
	integer userarg

	if (n1 .ge. n2) return

	nst=0
	k=n1
	l=n2

C sort sequence K...L

1	i=k
	j=l
	m=(k+l)/2
	call ld(m,userarg)			! load mid. element

C increment I until an element .GE. mid. element found

10	if (i.gt.l) goto 91
	call comp(i, tst,userarg)
	if (tst .ge. 0) goto 20
	i=i+1
	goto 10
91	stop 'QUICK_SORT: mistake in subroutine COMP, EXCH or LD'

C decrement J until an element .LT. mid. element found

20	if (j.lt.k) goto 91
	call comp(j, tst,userarg)
	if (tst .le. 0) goto 30
	j=j-1
	goto 20

30	if (i .ge. j) goto 40
	call exch(i, j,userarg)
	i=i+1
	j=j-1
	if (i .le. j) goto 10

C now all elements between K and I are .LE. mid. element
C and all elements between J and L are .GE. mid. element

40	i=i-1
	j=j+1
	if (i-k .gt. l-j) then		! K...I sequence is longer

	  if (i .gt. k) then		! if sequence longer than 1
	    nst=nst+1			! then push K,J to stack
	    if (nst .gt. nstack) goto 92
	    kst(nst)=k
	    lst(nst)=i
	  end if
	  k=j				! sort sequence J...L immediately

	else				! sequence J...L is longer

	  if (l .gt. j) then		! if sequence longer than 1
	    nst=nst+1			! then push J,L to stack
	    if (nst .gt. nstack) goto 92
	    kst(nst)=j
	    lst(nst)=l
	  end if
	  l=i				! sort sequence K...I immediately

	end if

	if (k .lt. l) goto 1		! if sequence longer than one then sort

90	if (nst .eq. 0) return		! no more sequences to sort
	k=kst(nst)			! pop K,L from stack
	l=lst(nst)
	nst=nst-1
	goto 1

92	stop 'QUICK_SORT: stack overflow'
	end