      function lleadz(i)
      implicit integer (a-z)
      parameter ( numbits=32 )
      logical bit

      lz=0
      do 1 l=1,numbits
      bit=btest(i,numbits-l)
      if(bit) then
        go to 2
      else
        lz=lz+1
      end if
    1 continue
    2 continue
      lleadz=lz

      return
      end

      function mmask(i)
      implicit integer (a-z)
      parameter ( numbits=64 )

      im=0

      if((i.gt.0).and.(i.le.63)) then

      do 1 l=1,i
      im=ibset(im,numbits-l)
    1 continue

      else if((i.ge.64).and.(i.le.128)) then

      do 2 l=1,128-i
      im=ibset(im,l-1)
    2 continue

      else

      im=0

      end if
      mmask=im
      return
      end

      subroutine ibatcher(xone,index,n)
      implicit integer (a-z)
      integer   xone
      dimension xone(n)
      dimension index(n)
      integer mmask,lleadz

      if(n.eq.1) return
C      P0=AND(N-1,MASK(LEADZ(N-1)+1))
C THIS PART TAKES CARE OF THE PREVIOUS LINE ***** FOR NON CRAY'S
      l=0
      nn=n-1
      do 2 m=1,64
      nn=nn/2
      if(nn.eq.0) go to 3
      l=l+1
    2 continue
    3 continue
      if(n.gt.1) then
        p0=2**l
      else
        p0=0
      end if
c--------------------------------
      p=p0
   20 continue
      q=p0
      r=0
      d=p
      i=0
   30 continue
      if(and(i,p).eq.r) then

cDIR$ IVDEP
      do 32 j=1,min(p,n-i-d)
      if(xone(index(j+i)).gt.xone(index(j+i+d))) then
        isave=index(j+i+d)
        index(j+i+d)=index(j+i)
        index(j+i)=isave
      end if
   32 continue

      end if

      i=i+p
      if((i+d).lt.n) go to 30
   50 continue
      if(q.ne.p) then
      d=q-p
      q=q/2
      r=p
      i=0
      go to 30
      end if
   60 continue
      p=p/2
      if(p.gt.4) go to 20
   70 continue
      q=p0
      r=0
      d=p
   80 continue
      do 82 i=1,p
cDIR$ IVDEP
      do 81 j=i,n-d-r,2*p
      if(xone(index(j+r)).gt.xone(index(j+r+d))) then
        isave=index(j+r+d)
        index(j+r+d)=index(j+r)
        index(j+r)=isave
      end if
   81 continue
   82 continue

      if(q.ne.p) then
      d=q-p
      q=q/2
      r=p
      go to 80
      end if

      p=p/2
      if(p.ne.0) go to 70

      return
      end

      subroutine batcher(xone,index,n)
      implicit integer (a-z)
      real*8   xone
      dimension xone(n)
      dimension index(n)
      integer mmask,lleadz

      if(n.eq.1) return
C      P0=AND(N-1,MASK(LEADZ(N-1)+1))
C THIS PART TAKES CARE OF THE PREVIOUS LINE ***** FOR NON CRAY'S
      l=0
      nn=n-1
      do 2 m=1,64
      nn=nn/2
      if(nn.eq.0) go to 3
      l=l+1
    2 continue
    3 continue
      if(n.gt.1) then
        p0=2**l
      else
        p0=0
      end if
c--------------------------------
      p=p0
   20 continue
      q=p0
      r=0
      d=p
      i=0
   30 continue
      if(and(i,p).eq.r) then

cDIR$ IVDEP
      do 32 j=1,min(p,n-i-d)
      if(xone(index(j+i)).gt.xone(index(j+i+d))) then
        isave=index(j+i+d)
        index(j+i+d)=index(j+i)
        index(j+i)=isave
      end if
   32 continue

      end if

      i=i+p
      if((i+d).lt.n) go to 30
   50 continue
      if(q.ne.p) then
      d=q-p
      q=q/2
      r=p
      i=0
      go to 30
      end if
   60 continue
      p=p/2
      if(p.gt.4) go to 20
   70 continue
      q=p0
      r=0
      d=p
   80 continue
      do 82 i=1,p
cDIR$ IVDEP
      do 81 j=i,n-d-r,2*p
      if(xone(index(j+r)).gt.xone(index(j+r+d))) then
        isave=index(j+r+d)
        index(j+r+d)=index(j+r)
        index(j+r)=isave
      end if
   81 continue
   82 continue

      if(q.ne.p) then
      d=q-p
      q=q/2
      r=p
      go to 80
      end if

      p=p/2
      if(p.ne.0) go to 70

      return
      end
