      subroutine heap_per(n,x,iper) ! bernd berg october 5 2000.  
c as heapsort, in addition the array  iper  keeps track of the permutation.
c permuted data are related to original by xout(i)=xin(iper(i)).
      include '../../libs/fortran/implicit.sta'
      dimension x(n),iper(n)
c
      do i=1,n
        iper(i)=i ! initialize iper.
      end do
c
c building the heap:
c ==================
      nhalf=n/2
      do m=nhalf,1,-1
        temp=x(m) 
        itmp=iper(m)
        i1=m
        i2=m+m
1       continue
c         index the larger of x(i2) and x(i2+1) as i2:
          if(i2.lt.n .and. x(i2+1).gt.x(i2)) i2=i2+1
          if(x(i2).gt.temp) then 
            x(i1)=x(i2)
            iper(i1)=iper(i2)
            i1=i2
            i2=i1+i1
          else  
            i2=n+1
          endif
        if(i2.le.n) go to 1
        x(i1)=temp ! final position for temp=x(m).
        iper(i1)=itmp
      end do
c
c sorting the heap:
c =================
      do i=n,3,-1 ! i is the number of still competing elements.
        temp=x(i)
        itmp=iper(i)
        x(i)=x(1) ! store top of the heap.
        iper(i)=iper(1)
        i1=1
        i2=2
2       continue
c         index the larger of x(i2) and x(i2+1) as i2:
          if((i2+1).lt.i .and. x(i2+1).gt.x(i2)) i2=i2+1
          if(x(i2).gt.temp) then  
            x(i1)=x(i2)
            iper(i1)=iper(i2)
            i1=i2
            i2=i1+i1
          else 
            i2=i
          endif 
        if(i2.lt.i) go to 2
        x(i1)=temp ! final position for temp=x(i).
        iper(i1)=itmp
      end do 
      temp=x(2)
      itmp=iper(2)
      x(2)=x(1)
      iper(2)=iper(1)
      x(1)=temp
      iper(1)=itmp
      return
      end
