website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/rank.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   box assisted sorting/ranking utilities
  23 c   author T. Schreiber (1998) based on earlier versions
  24 c===========================================================================
  25       subroutine rank(nmax,x,list)
  26 c  rank points in x
  27       parameter(nptr=100000)
  28       dimension x(nmax), list(nmax), jptr(0:nptr)
  29
  30       call minmax(nmax,x,xmin,xmax)
  31       if(xmin.eq.xmax) then
  32          do 10 n=1,nmax
  33  10         list(n)=n
  34          return
  35       endif
  36       nl=min(nptr,nmax/2)
  37       sc=(nl-1)/(xmax-xmin)
  38       do 20 i=0,nl
  39  20      jptr(i)=0
  40       do 30 n=1,nmax
  41          xn=x(n)
  42          i=int((xn-xmin)*sc)
  43          ip=jptr(i)
  44          if ((ip.eq.0).or.(xn.le.x(ip))) then
  45             jptr(i)=n
  46          else
  47  1          ipp=ip
  48             ip=list(ip)
  49             if ((ip.gt.0).and.(xn.gt.x(ip))) goto 1
  50             list(ipp)=n
  51          endif
  52  30      list(n)=ip
  53       n=0
  54       do 40 i=0,nl
  55          ip=jptr(i)
  56  2       if (ip.eq.0) goto 40
  57          n=n+1
  58          ipp=ip
  59          ip=list(ip)
  60          list(ipp)=n
  61          goto 2
  62 40       continue
  63       end
  64
  65       subroutine indexx(nmax,x,list)
  66 c make index table using rank
  67       dimension x(nmax), list(nmax)
  68
  69       call rank(nmax,x,list)
  70       call rank2index(nmax,list)
  71       end
  72
  73       subroutine rank2index(nmax,list)
  74 c converts a list of ranks into an index table (or vice versa) in place
  75       integer list(nmax)
  76
  77       do 10 n=1,nmax
  78  10      list(n)=-list(n)
  79       do 20 n=1,nmax
  80          if(list(n).gt.0) goto 20               ! has been put in place already
  81          ib=n
  82          im=-list(n)
  83  1       it=-list(im)
  84          list(im)=ib
  85          if(it.ne.n) then
  86             ib=im
  87             im=it
  88             goto 1
  89          else
  90             list(n)=im
  91          endif
  92  20      continue
  93       end
  94
  95       subroutine sort(nmax,x,list)
  96 c sort using rank and rank2sort
  97       dimension x(nmax), list(nmax)
  98
  99       call rank(nmax,x,list)
 100       call rank2sort(nmax,x,list)
 101       end
 102
 103       subroutine rank2sort(nmax,x,list)
 104 c sort x using list of ranks
 105       dimension x(nmax), list(nmax)
 106
 107       do 10 n=1,nmax
 108  10      list(n)=-list(n)
 109       do 20 n=1,nmax
 110          if(list(n).gt.0) goto 20               ! has been put in place already
 111          ib=n
 112          hb=x(n)
 113  1       it=-list(ib)
 114          list(ib)=it
 115          ht=x(it)
 116          x(it)=hb
 117          if(it.ne.n) then
 118             ib=it
 119             hb=ht
 120             goto 1
 121          endif
 122  20      continue
 123       end
 124
 125       subroutine index2sort(nmax,x,list)
 126 c sort x using list of indices
 127       dimension x(nmax), list(nmax)
 128
 129       do 10 n=1,nmax
 130  10      list(n)=-list(n)
 131       do 20 n=1,nmax
 132          if(list(n).gt.0) goto 20               ! has been put in place already
 133          ib=n
 134          h=x(n)
 135  1       it=-list(ib)
 136          list(ib)=it
 137          if(it.ne.n) then
 138             x(ib)=x(it)
 139             ib=it
 140             goto 1
 141          else
 142             x(ib)=h
 143          endif
 144  20      continue
 145       end
 146
 147       function which(nmax,x,k,list)
 148       dimension x(nmax), list(nmax)
 149
 150       call indexx(nmax,x,list)
 151       which=x(list(k))
 152       end
 153