ranl_skip2 

Routine to skip forward 2**d seeds for the RANL and RANL_NORMAL random number generators

FORMAT

  RANL_SKIP2 (d, s1, s2, s1_new, s2_new)

Arguments

  d                   integer*4
                      On entry, an integer d>0 specifying the number 2**d of
                      seeds to skip starting at (s1,s2) in the L'Ecuyer
                      algorithm.

  s1                  integer*4
                      On entry, a starting seed s1>=1.

  s2                  integer*4
                      On entry, a starting seed s2>=1.

  s1_new              integer*4
                      On exit, a new starting seed 2**d iterations away from
                      s1.

  s2_new              integer*4
                      On exit, a new starting seed 2**d iterations away from
                      s2.

Description

  The RANL_SKIP2 routine is used to get starting seeds for parallel,
  independent streams of random numbers. The new starting seeds are computed
  using the following algorithms:

  s1_new = a1**(2**d)*s1  mod m1
  s2_new = a2**(2**d)*s2  mod m2

  where a1,a2,m1,m2 are the constants defining the L'Ecuyer method.

  The following code example calls RANL_SKIP2 to set up separate seeds for 4
  streams, (nprocs=4), from RANL. It computes the averages per stream.

Example

        integer nprocs,n,hop
        parameter (nprocs=4)
        parameter (n=16,hop=14)
        integer*4 j,k
        real*4 v(n,nprocs)
        integer*4 s1val(nprocs),s2val(nprocs)
        real*4 sum1(nprocs)

  c     get seeds (2**hop apart) for separate streams
        s1val(1)=1
        s2val(1)=1
        do j=2,nprocs
        call ranl_skip2(hop,s1val(j-1),s2val(j-1),s1val(j),s2val(j))
        end do

  c     parallel calls to ranl can be done with KAP directives:
  C*$*  ASSERT CONCURRENT CALL
  C*$*  ASSERT DO (CONCURRENT)
        do j=1,nprocs
           call ranl(s1val(j),s2val(j),v(1,j),n)
           sum1(j)=0.0
           do k=1,n
              sum1(j)=sum1(j)+v(k,j)
           end do
        end do

        print*,' per-stream averages'
        do j=1,nprocs
           print*,sum1(j)/n
        end do

        end