c**        program temp
c**        implicit real*8 (a-h,o-z)
c**        real*8 x(10,2), nx(10,2)
c**        integer ia(10)
c**  c
c**  c
c**        n = 5 
c**        x(1,1) = 2.d0
c**        x(2,1) = 6.d0
c**        x(3,1) = 3.d0
c**        x(4,1) = 1.d0
c**        x(5,1) = 9.d0
c**  c
c**  c     have to assign 1..n to a, which is going to be the permutation
c**  c     order after sort.  Also have to set up a nx if you want to keep
c**  c     x unchanged.
c**  c
c**  c
c**        do 10 i = 1,n  
c**         x(i,2) = x(6-i,1)
c**         nx(i,1) = x(i,1)
c**         nx(i,2) = x(i,2)
c**         ia(i) = i
c**         print *, (x(i,j),j=1,2)
c**   10   continue
c**  c
c**  c
c**  c     will sort the 1st column and obtain a which can be used for other
c**  c     columns
c**  c
c**  c
c**        call sort(nx(1,1),ia,1,5,5)
c**  c
c**        do 20 i = 1,n
c**   20   nx(i,2) = x(ia(i),2)
c**  c
c**        print *, (ia(i),i=1,5)
c**  c
c**  c
c**        do 30 i = 1,n
c**   30    print *, (nx(i,j),j=1,2)
c**  c
c**        stop
c**        end
      

c*******************************************************************
      subroutine sort (v,ia,ii,jj,n)
c
c     puts into a the permutation vector which sorts v into
c     increasing order.  only elements from ii to jj are considered.
c     arrays iu(k) and il(k) permit sorting up to 2**(k+1)-1 elements
c
c     this is a modification of cacm algorithm #347 by r. c. singleton,
c     which is a modified hoare quicksort.
c
      implicit real*8 (a-h,o-z)
      real*8  v(n),iu(5000),il(5000)
      integer mmbb,jjuu
      integer ia(n)
      m=1
      i=ii
      j=jj
 10   if (i.ge.j) go to 80
 20   k=i
      ij=(j+i)/2
      mmbb=ia(ij)
      vss=v(ij)
      if (v(i).le.vss) go to 30
      ia(ij)=ia(i)
      ia(i)=mmbb
      mmbb=ia(ij)
      v(ij)=v(i)
      v(i)=vss
      vss=v(ij)
 30   l=j
      if (v(j).ge.vss) go to 50
      ia(ij)=ia(j)
      ia(j)=mmbb
      mmbb=ia(ij)
      v(ij)=v(j)
      v(j)=vss
      vss=v(ij)
      if (v(i).le.vss) go to 50
      ia(ij)=ia(i)
      ia(i)=mmbb
      mmbb=ia(ij)
      v(ij)=v(i)
      v(i)=vss
      vss=v(ij)
      go to 50
 40   ia(l)=ia(k)
      ia(k)=jjuu
      v(l)=v(k)
      v(k)=vrr
 50   l=l-1
      if (v(l).gt.vss) go to 50
      jjuu=ia(l)
      vrr=v(l)
 60   k=k+1
      if (v(k).lt.vss) go to 60
      if (k.le.l) go to 40
      if (l-i.le.j-k) go to 70
      il(m)=i
      iu(m)=l
      i=k
      m=m+1
      go to 90
 70   il(m)=k
      iu(m)=j
      j=l
      m=m+1
      go to 90
 80   m=m-1
      if (m.eq.0) return
      i=il(m)
      j=iu(m)
 90   if (j-i.gt.10) go to 20
      if (i.eq.ii) go to 10
      i=i-1
 100  i=i+1
      if (i.eq.j) go to 80
      mmbb=ia(i+1)
      vss=v(i+1)
      if (v(i).le.vss) go to 100
      k=i
 110  ia(k+1)=ia(k)
      v(k+1)=v(k)
      k=k-1
      if (vss.lt.v(k)) go to 110
      ia(k+1)=mmbb
      v(k+1)=vss
      go to 100
      end