indexxd.for

Go to the documentation of this file.

c> \ingroup core
c> Subroutine to perform ?? indexing on floating point data (double precision)
c> 
c> \param[in] n number of iterations (rows?)
c> \param[in] nd  array and index dimensions
c> \param[in] arr input data array
c> \param[out] indx index out
c>
c> ## Changelog
c>
c> ### 1/8/2004: 
c> Modified to allow `indx` and `arr` to
c> be DIMENSIONED differently than the number of good
c> values they contain
c>
c> ### 11/7/2003
c> Modified to REAL*8 by D. Smith, 
c>
      SUBROUTINE indexxd(n,nd,arr,indx)
c - Modified to REAL*8 by D. Smith, 11/7/2003
c - Further modified 1/8/2004 to allow indx and arr to
c - be DIMENSIONED differently than the number of good
c - values they contain

      integer nd
      INTEGER n,indx(nd),M,NSTACK
      REAL*8 arr(nd)
      parameter(m=7,nstack=50)
      INTEGER i,indxt,ir,itemp,j,jstack,k,l,istack(nstack)
      REAL*8 a
      do 11 j=1,n
        indx(j)=j
11    continue
      jstack=0
      l=1
      ir=n
1     if(ir-l.lt.m)then
        do 13 j=l+1,ir
          indxt=indx(j)
          a=arr(indxt)
          do 12 i=j-1,1,-1
            if(arr(indx(i)).le.a)goto 2
            indx(i+1)=indx(i)
12        continue
          i=0
2         indx(i+1)=indxt
13      continue
        if(jstack.eq.0)return
        ir=istack(jstack)
        l=istack(jstack-1)
        jstack=jstack-2
      else
        k=(l+ir)/2
        itemp=indx(k)
        indx(k)=indx(l+1)
        indx(l+1)=itemp
        if(arr(indx(l+1)).gt.arr(indx(ir)))then
          itemp=indx(l+1)
          indx(l+1)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l)).gt.arr(indx(ir)))then
          itemp=indx(l)
          indx(l)=indx(ir)
          indx(ir)=itemp
        endif
        if(arr(indx(l+1)).gt.arr(indx(l)))then
          itemp=indx(l+1)
          indx(l+1)=indx(l)
          indx(l)=itemp
        endif
        i=l+1
        j=ir
        indxt=indx(l)
        a=arr(indxt)
3       continue
          i=i+1
        if(arr(indx(i)).lt.a)goto 3
4       continue
          j=j-1
        if(arr(indx(j)).gt.a)goto 4
        if(j.lt.i)goto 5
        itemp=indx(i)
        indx(i)=indx(j)
        indx(j)=itemp
        goto 3
5       indx(l)=indx(j)
        indx(j)=indxt
        jstack=jstack+2
        if(jstack.gt.nstack)pause 'NSTACK too small in indexx'
        if(ir-i+1.ge.j-l)then
          istack(jstack)=ir
          istack(jstack-1)=i
          ir=j-1
        else
          istack(jstack)=j-1
          istack(jstack-1)=l
          l=i
        endif
      endif
      goto 1
      END