c
c  CVS Info
c  $Date: 2005/01/10 21:56:21 $
c  $Revision: 1.2 $
c  $RCSfile: newrandom.f,v $
c  $Name: rel_5 $
c
      subroutine initiran(seqnum,gseed,state)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine initializes the state for the seqnum(th) member        c
c  of the family of pseudorandom sequences.  initiran must             c
c  be called once (and only once) for each sequence from which         c
c  random binary or integer values are wanted.                         c
c                                                                      c
c  l, k, s, and mbits are properties associated with the family        c
c  of generators used by this package.  l and k represent a            c
c  particular primitive trinomial (mod 2), s is the "canonical bit"    c
c  associated with that trinomial, and mbits is the number of          c
c  mantissa bits used by the floating point pseudorandom generator.    c
c  for l = 158, state must have a declared dimension of at least       c
c  4 x 158 = 632 in the calling program.                               c
c                                                                      c
c  After initializing a sequence with initiran, multiple calls         c
c  can be made to iranx or brand64,  and the state vector will         c
c  stay current.   The state vector must not be modified by the user   c
c  code.   If more than one sequence is to be active at any time,      c
c  there must be a searate state vector for each active sequence.      c
c                                                                      c
c  gseed is the integer value of the user-supplied "global seed."      c
c  The global seed is intended to be fixed for the duration of the     c
c  run and should be the same value for every sequence on every        c
c  processor.   It can be used to modify the random data, in a         c
c  predictable way, on a per-run basis.  seqnum must be in the         c
c  range {0,1,2,3,..., 65535 = 2^16 - 1 = z'0000ffff'} and gseed       c
c  must be in the range {0,1,2,3,..., 65534 = 2^16 - 2 = z'0000fffe'}. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer seqnum,gseed,state(4*l)
c----------------------------------------------------------------------c
      integer hi,lo,a,q,r,t31m1,seed,i
      integer treg1,treg2,highbit1,highbit2
      parameter(a = 16807, q = 127773, r = 2836, t31m1 = z'7fffffff')
      
C--CVS variable declaration
      TYPE CVS
      sequence
      character( 160 )     string
      integer              stringend
      END TYPE CVS       
C--CVS initilaize variables    
      TYPE( CVS ),save :: CVS_INFO =
     $  CVS("BMARKGRP $Date: 2005/01/10 21:56:21 $ $Revision: 1.2 $" //
     $      "$RCSfile: newrandom.f,v $ $Name: rel_5 $", 0)
c
      call tausinit(seqnum,gseed,treg1,treg2)
c
c  initialize first state register:
      seed     = iand(treg1,t31m1)
      highbit1 = ieor(seed,treg1)
      state(l) = 0
      do  i=1,l-1
         hi        = seed/q
         lo        = seed - (hi*q)
         seed      = a*lo - r*hi
         if(seed.le.0) seed = seed + t31m1
         state(i)  = seed + seed
      enddo
c
c  incorporate the high bit of the tuas. reg. in the
c  fibonacci register, in order to maintain uniqueness:
      if(highbit1.eq.0) then
         state(l-1) = iand(state(l-1),t31m1)
      else
         state(l-1) = ior (state(l-1),highbit1)
      endif
      state(s+1)   = ior(state(s+1),1)
c
c  initialize second state register:
      seed       = iand(treg2,t31m1)
      highbit2   = ieor(seed,treg2)
      state(2*l) = 0
      do  i=l+1,2*l-1
         hi        = seed/q
         lo        = seed - (hi*q)
         seed      = a*lo - r*hi
         if(seed.le.0) seed = seed + t31m1
         state(i)  = seed + seed
      enddo
c
c  incorporate the high bit of the tuas. reg. in the
c  fibonacci register, in order to maintain uniqueness:
      if(highbit2.eq.0) then
         state(2*l-1) = iand(state(2*l-1),t31m1)
      else
         state(2*l-1) = ior (state(2*l-1),highbit2)
      endif
      state(l+s+1) = ior(state(l+s+1),1)
      return
      end
c
c
c
c
c
      subroutine initfran(seqnum,gseed,state)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine initializes the state for the seqnum(th) member        c
c  of the family of pseudorandom sequences.  initfran must             c
c  be called once (and only once) for each sequence from which         c
c  random floating point values are wanted.                            c
c                                                                      c
c  l, k, s, and mbits are properties associated with the family        c
c  of generators used by this package.  l and k represent a            c
c  particular primitive trinomial (mod 2), s is the "canonical bit"    c
c  associated with that trinomial, and mbits is the number of          c
c  mantissa bits used by the floating point pseudorandom generator.    c
c  for l = 158, state must have a declared dimension of at least       c
c  2 x 158 = 316 in the calling program.                               c
c                                                                      c
c  After initializing a sequence with initfran, multiple calls         c
c  can be made to frand,  and the state vector will stay current.      c
c  The state vector must not be modified by the user code.   If        c
c  more than one sequence is to be active at any time, there must      c
c  be a searate state vector for each active sequence.                 c
c                                                                      c
c  gseed is the integer value of the user-supplied "global seed."      c
c  The global seed is intended to be fixed for the duration of the     c
c  run and should be the same value for every sequence on every        c
c  processor.   It can be used to modify the random data, in a         c
c  predictable way, on a per-run basis.  seqnum must be in the         c
c  range {0,1,2,3,..., 65535 = 2^16 - 1 = z'0000ffff'} and gseed       c
c  must be in the range {0,1,2,3,..., 65534 = 2^16 - 2 = z'0000fffe'}. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer seqnum,gseed
      real    state(2*l)
c----------------------------------------------------------------------c
      integer hi,lo,a,q,r,t31m1,seed,i,mask
      integer treg,tregx,highbit,t30m1
      parameter(a = 16807, q = 127773, r = 2836 , t31m1 = z'7fffffff')
      parameter(t30m1 = z'3fffffff')
      real    scale
c
      call tausinit(seqnum,gseed,treg,tregx)
      scale    = 1.
      do  i=1,mbits
         scale = scale * 0.5
      enddo
      mask     = ishft(1,mbits) - 2
      seed     = iand(treg,t31m1)
      highbit  = ieor(seed,treg)
c
c  initialize register:
      state(1) = scale*float(seqnum + seqnum)
      state(l) = 0.
      do  i=2,l-2
         hi    = seed/q
         lo    = seed - (hi*q)
         seed  = a*lo - r*hi
         if(seed.le.0) seed = seed + t31m1
         state(i) = scale*float(iand(ishft(seed,-(31-mbits)), mask))
      enddo
c
      hi    = seed/q
      lo    = seed - (hi*q)
      seed  = a*lo - r*hi
      if(seed.le.0) seed = seed + t31m1
      if(highbit.eq.0) then
         seed = iand(seed,t30m1)
      else
         seed = ior (seed,ishft(highbit,-1))
      endif
      state(l-1) = scale*float(iand(ishft(seed,-(31-mbits)), mask))
c
      state(s+1)  = state(s+1) + scale
      return
      end
c
c
c
c
c
      subroutine irand(n,state,x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine returns n 32-bit random integers in the array, x.      c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer n
      integer state(l,2,0:1),x(n)
c----------------------------------------------------------------------c
      integer i,old,new,xptr,kk,t31m1,t32m2
      parameter(t31m1 = z'7fffffff', t32m2 = z'fffffffe')
c
      old  = 1
      new  = 0
      xptr = 0
   10 if(xptr.ge.n) go to 99
         old  = new
         new  = ieor(new,1)
         kk = min(k,n-xptr)
         do  i=1,kk
            state(i,1,new) = state(k-kk+i,1,old) + state(l-kk+i,1,old)
            state(i,2,new) = state(k-kk+i,2,old) + state(l-kk+i,2,old)
         enddo
         do  i=kk+1,l
            state(i,1,new) = state(i-kk,1,old)
            state(i,2,new) = state(i-kk,2,old)
         enddo
         do  i=xptr+1, min(xptr+kk,n)
            x(i) = iand(ishft(state(kk+1 - (i-xptr),1,new),   -1),t31m1)
            x(i) = ieor( iand(state(kk+1 - (i-xptr),2,new),t32m2), x(i))
         enddo
         xptr = xptr + kk
         go to 10
   99 continue
c
c  update "old" state for subsequent calls:
      if(new.eq.1) then
         do  i=1,l
            state(i,1,old) = state(i,1,new)
            state(i,2,old) = state(i,2,new)
         enddo
      endif
      return
      end
c
c
c
c
c
      subroutine frand(n,state,x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine returns n random real numbers in the array, x.         c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer n
      real    state(l,0:1),x(n)
c----------------------------------------------------------------------c
      integer i,old,new,xptr,kk
c
      old  = 1
      new  = 0
      xptr = 0
   10 if(xptr.ge.n) go to 99
         old  = new
         new  = ieor(new,1)
         kk   = min(k,n-xptr)
         do  i=1,kk
            state(i,new) = state(k-kk+i,old) + state(l-kk+i,old)
         enddo
         do  i=1,kk
            if(state(i,new).ge.1.0) state(i,new) = state(i,new) - 1.0
         enddo
         do  i=kk+1,l
            state(i,new) = state(i-kk,old)
         enddo
         do  i=xptr+1, min(xptr+kk,n)
            x(i) = state(kk+1 - (i-xptr),new)
         enddo
         xptr = xptr + kk
         go to 10
   99 continue
c
c  update "old" state for subsequent calls:
      if(new.eq.1) then
         do  i=1,l
            state(i,old) = state(i,new)
         enddo
      endif
      return
      end
c
c
c
c
c
      subroutine brand64(n,state,x)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Used for producing random bit streams on machines with 64-bit       c
c  integers, this routine returns n words of random 64-bit integers    c
c  in the array, x.  This is the same string of bits that would        c
c  be produced on a 32-bit machine by calling irand with a             c
c  (2 x n)-long x vector.  Thus, for creating binary streams, use      c
c  irand on 32-bit machines, and brand64 on 64-bit machines.           c
c                                                                      c
c  The state vector is updated to continue producing random numbers    c
c  from the current sequence by subsequent calls to this routine.      c
c  The state vector should not be modified by the calling routine,     c
c  unless no more values are required from the current sequence.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer wrklen
c  wrklen must be even!!
      parameter(wrklen = 8192)
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer   n
      integer   state(4*l)
c
c  for some machines (e.g. SGI Power Challenge), the next line
c  should read "integer*8" instead of just "integer".
      integer*8 x(n),wrk8(wrklen),low32
c----------------------------------------------------------------------c
      integer   xptr,nleft,nleftd2,i
      integer   wrk(wrklen)
      data low32/ z'00000000FFFFFFFF'/
c
      xptr  = 1
      nleft = n
   10 if(nleft.le.0) go to 99
      nleftd2 = nleft/2
      if(2*nleft.gt.wrklen) then
         call irand(wrklen,state,wrk)
         do  i=1,wrklen
            wrk8(i) = wrk(i)
         enddo
         do  i=1,wrklen
            wrk8(i) = iand(wrk8(i),low32)
         enddo
         do  i=xptr,xptr+(wrklen/2)-1
            x(i) = ior(        wrk8(2*(i+1-xptr) - 1) ,
     $                  ishft( wrk8(2*(i+1-xptr)    ) , 32))
         enddo
         xptr  = xptr  + (wrklen/2)
         nleft = nleft - (wrklen/2)
      else
         call irand(2*nleft,state,wrk)
         do  i=1,2*nleft
            wrk8(i) = wrk(i)
         enddo
         do  i=1,2*nleft
            wrk8(i) = iand(wrk8(i),low32)
         enddo
         do  i=xptr,xptr+nleft-1
            x(i) = ior(        wrk8(2*(i+1-xptr) - 1) ,
     $                  ishft( wrk8(2*(i+1-xptr)    ) , 32))
         enddo
         xptr  = xptr + nleft
         nleft = 0
      endif
      go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine get_penum(seqnum,npes,nseqs,penum)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  given  a sequence number (seqnum), the number of pes (npes), and    c
c  the number of sequences (nseqs) this routine determines the number  c
c  of the pe on which sequence number seqnum is located.               c
c  Assumes that penum is in the range {0,1,2,...,npes-1} and that      c
c  sequences are numbered starting from 0.                             c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum,npes,nseqs,penum
c----------------------------------------------------------------------c
      integer q,r
c
      q     = nseqs/npes
      r     = nseqs - npes*q
      penum = seqnum/(q+1)
      if(seqnum.ge. ((q+1)*r) ) penum = r + (seqnum - ((q+1)*r))/q
      return
      end
c
c
c
c
c
      subroutine get_seq_nums(penum,npes,nseqs,seqfirst,seqlast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  This routine computes the indices of the first and last sequences   c
c  that will be associated with pe number penum, when there are npes   c
c  pes, and nseqs sequences.  Assumes that penum is in the range       c
c  {0,1,2,...,npes-1} and that sequences are numbered starting from 0. c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer penum,npes,nseqs,seqfirst,seqlast
c----------------------------------------------------------------------c
      integer q,r
c
      q        = nseqs/npes
      r        = nseqs - npes*q
      seqfirst = penum*q + min(penum,r)
      seqlast  = seqfirst + q - 1
      if(penum.lt.r) seqlast = seqlast + 1
      return
      end
c
c
c
c
c
      subroutine get_seq_bounds(seqnum,nseqs,lendata,chnksize,
     $                          datafirst,datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Given a particular sequence number (seqnum), the number of          c
c  sequences (nseqs), and the total length of the data (lendata) ,     c
c  this will compute datafirst and datalast, the positions in the      c
c  overall data set that sequence number seqnum occupies.   This       c
c  routine will force each sequence to have an integer number of data  c
c  values of length chnksize (with the exception of the last           c
c  sequence, which will have any partial chunk left over from a total  c
c  data length of lendata).  Assumes that data arrays are numbered     c
c  starting from 1, and that sequences are numbered starting from 0.   c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum,nseqs,lendata,datafirst,datalast,chnksize
c----------------------------------------------------------------------c
      integer q,r,nchunks
c
      nchunks   = lendata/chnksize
      q         = nchunks/nseqs
      r         = nchunks - nseqs*q
      datafirst = chnksize*(seqnum*q + min(seqnum,r)) + 1
      datalast  = datafirst + q*chnksize - 1
      if(seqnum.lt.r        ) datalast = datalast + chnksize
      if(seqnum.eq.(nseqs-1)) datalast = lendata
      return
      end
c
c
c
c
c
      subroutine get_my_idata(penum,npes,nseqs,lendata,chnksize,
     $                        gseed,x,datafirst,datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the integer data belonging to the PE numbered penum,        c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata           c
c  integers, a chunk size of chnksize integers, and a global seed      c
c  value of gseed.  The data belonging to PE penum will be returned in c
c  the integer array, x, starting with x(1),  but will be used by the  c
c  calling program as data items datafirst through datalast in the     c
c  overall lendata-long set that resides across all PEs.  For this     c
c  routine to work correctly, npes must be less than the total number  c
c  of chunks in the data set.  I.e. npes must be less than             c
c  lendata/chnksize.  The size of x must be declared to be at least as c
c  large as the value of maxlen returned by the routine check_maxlen.  c
c  Of course, this will be about lendata/npes, but the value of maxlen c
c  returned by check_maxlen is the exact value needed by the the       c
c  worst-case PE.                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer penum,npes,nseqs,lendata,chnksize,gseed,
     $        x(lendata),datafirst,datalast
c----------------------------------------------------------------------c
      integer state(4*l)
      integer nchnks,q,lendatac,r,ch_last,ch_frst,mynchnks,
     $        mydata_len,dlen,dlast,dfrst,seq_frst,nskip,ix,mydata_left,
     $        len_data1,nseq_len2,nseq_len1,seq_frst2,datafirst2,
     $        len_data2,seq_len1,seq_len2c,seq_len2,rc,seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      ch_last     = (penum+1)*q + min((penum+1),r)
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      mynchnks    = ch_last - ch_frst + 1
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if(penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if(datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
      call initiran(seq_frst,gseed,state)
      call irand(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
      call irand(dlen,state,x(ix))
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if(mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast - dfrst + 1,mydata_left)
         call initiran(seq_frst,gseed,state)
         call irand(dlen,state,x(ix))
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine get_my_fdata(penum,npes,nseqs,lendata,chnksize,
     $                        gseed,x,datafirst,datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the floating point data belonging to the PE numbered penum, c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata           c
c  integers, a chunk size of chnksize reals, and a global seed         c
c  value of gseed.  The data belonging to PE penum will be returned in c
c  the real array, x, starting with x(1),  but will be used by the     c
c  calling program as data items datafirst through datalast in the     c
c  overall lendata-long set that resides across all PEs.  For this     c
c  routine to work correctly, npes must be less than the total number  c
c  of chunks in the data set.  I.e. npes must be less than             c
c  lendata/chnksize.  The size of x must be declared to be at least as c
c  large as the value of maxlen returned by the routine check_maxlen.  c
c  Of course, this will be about lendata/npes, but the value of maxlen c
c  returned by check_maxlen is the exact value needed by the the       c
c  worst-case PE.                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer penum,npes,nseqs,lendata,chnksize,gseed,
     $        datafirst,datalast
      real    x(lendata)
c----------------------------------------------------------------------c
      real    state(2*l)
      integer nchnks,q,lendatac,r,ch_last,ch_frst,mynchnks,
     $        mydata_len,dlen,dlast,dfrst,seq_frst,nskip,ix,mydata_left,
     $        len_data1,nseq_len2,nseq_len1,seq_frst2,datafirst2,
     $        len_data2,seq_len1,seq_len2c,seq_len2,rc,seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      ch_last     = (penum+1)*q + min((penum+1),r)
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      mynchnks    = ch_last - ch_frst + 1
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if(penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if(datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
      call initfran(seq_frst,gseed,state)
      call frand(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
      call frand(dlen,state,x(ix))
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if(mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast - dfrst + 1,mydata_left)
         call initfran(seq_frst,gseed,state)
         call frand(dlen,state,x(ix))
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine get_my_bdata64(penum,npes,nseqs,lendata,chnksize,
     $                        gseed,x,datafirst,datalast)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  returns the 64-bit integer data belonging to the PE numbered penum, c
c  given that there are npes PEs and given that the user has           c
c  specified nseqs sequences, a total data length of lendata 64-bit    c
c  integers, a chunk size of chnksize 64-bit integers, and a global    c
c  seed value of gseed.  The data belonging to PE penum will be        c
c  returned in the 64-bit integer array, x, starting with x(1),  but   c
c  will be used by the calling program as data items datafirst         c
c  through datalast in the overall lendata-long set that resides       c
c  across all PEs.  For this routine to work correctly, npes must be   c
c  less than the total number of chunks in the data set.  I.e. npes    c
c  must be less than lendata/chnksize.  The size of x must be          c
c  declared to be at least as large as the value of maxlen returned    c
c  by the routine check_maxlen.  Of course, this will be about         c
c  lendata/npes, but the value of maxlen returned by check_maxlen is   c
c  the exact value needed by the the worst-case PE.                    c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer l,k,s,mbits
      parameter(l = 158, k = 128, s =  63, mbits = 23)
      integer   penum,npes,nseqs,lendata,chnksize,gseed,
     $          datafirst,datalast
c
c  for some machines (e.g. SGI Power Challenge), the next line
c  should read "integer*8" instead of just "integer".
      integer*8 x(lendata)
c----------------------------------------------------------------------c
      integer state(4*l)
      integer nchnks,q,lendatac,r,ch_last,ch_frst,mynchnks,
     $        mydata_len,dlen,dlast,dfrst,seq_frst,nskip,ix,mydata_left,
     $        len_data1,nseq_len2,nseq_len1,seq_frst2,datafirst2,
     $        len_data2,seq_len1,seq_len2c,seq_len2,rc,seq_len1c
c
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      ch_last     = (penum+1)*q + min((penum+1),r)
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      mynchnks    = ch_last - ch_frst + 1
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if(penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
c
      lendatac = lendata/chnksize
      seq_len2c = lendatac/nseqs
      seq_len1c = seq_len2c + 1
      rc        = lendatac - seq_len2c*nseqs
c
      seq_len2  = chnksize*seq_len2c
      seq_len1  = chnksize*seq_len1c
      nseq_len1 = rc
      nseq_len2 = nseqs - rc
c
      len_data1 = nseq_len1*seq_len1
      len_data2 = nseq_len2*seq_len2
c
      if(datafirst.gt.len_data1) then
         datafirst2 = datafirst - len_data1
         seq_frst2  = datafirst2/seq_len2
         nskip      = datafirst2 - seq_frst2*seq_len2 - 1
         seq_frst   = seq_frst2 + nseq_len1
      else
         seq_frst   = datafirst/seq_len1
         nskip      = datafirst - seq_frst*seq_len1 - 1
      endif
c
      call get_seq_bounds(seq_frst,nseqs,lendata,chnksize,dfrst,dlast)
      dlen        = dlast - dfrst + 1
      call initiran(seq_frst,gseed,state)
      call brand64(nskip,state,x)
c
      dlen        = dlen - nskip
      mydata_left = mydata_len
      ix          = 1
      call brand64(dlen,state,x(ix))
      mydata_left = mydata_left - dlen
      ix          = ix + dlen
c
   10 if(mydata_left.le.0) go to 99
         seq_frst = seq_frst + 1
         call get_seq_bounds(seq_frst,nseqs,lendata,
     $                       chnksize,dfrst,dlast)
         dlen     = min(dlast - dfrst + 1,mydata_left)
         call initiran(seq_frst,gseed,state)
         call brand64(dlen,state,x(ix))
         mydata_left = mydata_left - dlen
         ix          = ix + dlen
         go to 10
   99 continue
      return
      end
c
c
c
c
c
      subroutine check_maxlen(npes,nseqs,lendata,chnksize,maxlen)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer npes,nseqs,lendata,chnksize,maxlen
c----------------------------------------------------------------------c
      integer mynchnks,ch_frst,datafirst,mydata_len,datalast,ch_last,
     $        nchnks,penum,q,r,lendatac
c
      penum       = 0
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      ch_last     = (penum+1)*q + min((penum+1),r)
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      mynchnks    = ch_last - ch_frst + 1
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if(penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
      maxlen      = mydata_len
c
      penum       = npes-1
      nchnks      = lendata/chnksize
      q           = nchnks/npes
      lendatac    = nchnks*chnksize
      r           = nchnks - npes*q
      ch_last     = (penum+1)*q + min((penum+1),r)
      ch_frst     =  penum   *q + min( penum   ,r) + 1
      mynchnks    = ch_last - ch_frst + 1
      datafirst   = chnksize*(ch_frst-1) + 1
      datalast    = chnksize*(ch_last-1) + chnksize
      if(penum.eq.npes-1) datalast   = lendata
      mydata_len  = datalast - datafirst + 1
      maxlen      = max(maxlen,mydata_len)
c
      return
      end
c
c
c
c
c
      subroutine tausinit(seqnum,gseed,treg64,treg96)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  uses a 32-bit tausworthe generator to map the pair (seqnum,gseed)   c
c  into a 32-bit pseudorandom number, for subsequent use in filling    c
c  out the fibonacci register.  This routine produces two values,      c
c  treg64 and treg96, which represent 64 and 96 steps of the           c
c  initial fill.  Since the calling routine will only be using         c
c  the low 31 bits (mod 2^31 - 1), there are 6 cases to avoid.         c
c  The bad values for the tausworthe register are (in order):          c
c    T^(-64)(80000000),T^(-64)(7fffffff),T^(-64)(ffffffff),            c
c    T^(-96)(80000000),T^(-96)(7fffffff),T^(-96)(ffffffff).            c
c  these cases correspond to (gseed,seqnum) pairs (in decimal):        c
c  (43553,65535), (26141,65535), (52283,65535),                        c
c  (50425,37887), (17320,35839), (34642, 6143).                        c
c  these are mapped to multiples of the prime number 9349, which       c
c  is approximately 65536/7 and which do not correspond to valid       c
c  initial values of seqnum and gseed that would be passed into        c
c  this routine.                                                       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer seqnum,gseed,treg64,treg96
c----------------------------------------------------------------------c
      integer seqmask,gshft,i,treg
      parameter(seqmask = z'0000ffff')
      parameter(gshft   = 16)
      integer tbad(6),t_ok(6)
      data tbad/ z'AA220000',z'661E0000',z'CC3C0000',
     $           z'C4FA9400',z'43A98C00',z'87531800'/
      data t_ok/ z'00002485',z'0000490a',z'00006d8f',
     $           z'00009214',z'0000b699',z'0000db1e'/
c
      treg     = ishft(gseed + 1,gshft)
      treg     = ior(treg,iand(seqnum + 1,seqmask))
      do  i=1,6
         if(treg.eq.tbad(i)) treg = t_ok(i)
      enddo
c
      call taus32(treg,64)
      treg64 = treg
      call taus32(treg,32)
      treg96 = treg
      return
      end
c
c
c
c
c
      subroutine taus32(treg,nsteps)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  steps the 32-bit tausworthe generator nsteps times, using the       c
c  mod 2 primitive polynomial (32,7,5,3,2,1,0).  The table entries     c
c  are used to look up the parity of the least significant 7 bits.     c
c  Tausworthe generators are discussed in Numerical Recipes by         c
c  Press, et al.  Cambridge University Press, c1986, pp 209-213.       c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      implicit none
      integer treg,nsteps
c----------------------------------------------------------------------c
      integer i,t31,t32m1,table0(0:127),table1(0:127),bit31,low7
      parameter(t31     = z'80000000')
      parameter(t32m1   = z'FFFFFFFF')
c
      data table0 /
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1/
      data table1 /
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0, 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0,
     $ 0,1,1,0,1,0,0,1,0,1,1,0,1,0,0,1, 1,0,0,1,0,1,1,0,1,0,0,1,0,1,1,0/
c
      do  i=1,nsteps
         bit31 = iand(treg,t31)
         low7  = iand(treg,127)
         treg  = iand(ishft(treg,1),t32m1)
         if(bit31.eq.0) then
            treg = ior(treg,table0(low7))
         else
            treg = ior(treg,table1(low7))
         endif
      enddo
      return
      end