examples/smpi/NAS/common/randdpvec.f

   1 c---------------------------------------------------------------------
   2       double precision function randlc (x, a)
   3 c---------------------------------------------------------------------
   4
   5 c---------------------------------------------------------------------
   6 c
   7 c   This routine returns a uniform pseudorandom double precision number in the
   8 c   range (0, 1) by using the linear congruential generator
   9 c
  10 c   x_{k+1} = a x_k  (mod 2^46)
  11 c
  12 c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  13 c   before repeating.  The argument A is the same as 'a' in the above formula,
  14 c   and X is the same as x_0.  A and X must be odd double precision integers
  15 c   in the range (1, 2^46).  The returned value RANDLC is normalized to be
  16 c   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
  17 c   the new seed x_1, so that subsequent calls to RANDLC using the same
  18 c   arguments will generate a continuous sequence.
  19 c
  20 c   This routine should produce the same results on any computer with at least
  21 c   48 mantissa bits in double precision floating point data.  On 64 bit
  22 c   systems, double precision should be disabled.
  23 c
  24 c   David H. Bailey     October 26, 1990
  25 c
  26 c---------------------------------------------------------------------
  27
  28       implicit none
  29
  30       double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
  31       parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
  32      >  t46 = t23 ** 2)
  33
  34 c---------------------------------------------------------------------
  35 c   Break A into two parts such that A = 2^23 * A1 + A2.
  36 c---------------------------------------------------------------------
  37       t1 = r23 * a
  38       a1 = int (t1)
  39       a2 = a - t23 * a1
  40
  41 c---------------------------------------------------------------------
  42 c   Break X into two parts such that X = 2^23 * X1 + X2, compute
  43 c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
  44 c   X = 2^23 * Z + A2 * X2  (mod 2^46).
  45 c---------------------------------------------------------------------
  46       t1 = r23 * x
  47       x1 = int (t1)
  48       x2 = x - t23 * x1
  49
  50
  51       t1 = a1 * x2 + a2 * x1
  52       t2 = int (r23 * t1)
  53       z = t1 - t23 * t2
  54       t3 = t23 * z + a2 * x2
  55       t4 = int (r46 * t3)
  56       x = t3 - t46 * t4
  57       randlc = r46 * x
  58       return
  59       end
  60
  61
  62
  63 c---------------------------------------------------------------------
  64 c---------------------------------------------------------------------
  65
  66       subroutine vranlc (n, x, a, y)
  67
  68 c---------------------------------------------------------------------
  69 c---------------------------------------------------------------------
  70
  71 c---------------------------------------------------------------------
  72 c   This routine generates N uniform pseudorandom double precision numbers in
  73 c   the range (0, 1) by using the linear congruential generator
  74 c
  75 c   x_{k+1} = a x_k  (mod 2^46)
  76 c
  77 c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  78 c   before repeating.  The argument A is the same as 'a' in the above formula,
  79 c   and X is the same as x_0.  A and X must be odd double precision integers
  80 c   in the range (1, 2^46).  The N results are placed in Y and are normalized
  81 c   to be between 0 and 1.  X is updated to contain the new seed, so that
  82 c   subsequent calls to RANDLC using the same arguments will generate a
  83 c   continuous sequence.
  84 c
  85 c   This routine generates the output sequence in batches of length NV, for
  86 c   convenience on vector computers.  This routine should produce the same
  87 c   results on any computer with at least 48 mantissa bits in double precision
  88 c   floating point data.  On Cray systems, double precision should be disabled.
  89 c
  90 c   David H. Bailey    August 30, 1990
  91 c---------------------------------------------------------------------
  92
  93       integer n
  94       double precision x, a, y(*)
  95
  96       double precision r23, r46, t23, t46
  97       integer nv
  98       parameter (r23 = 2.d0 ** (-23), r46 = r23 * r23, t23 = 2.d0 ** 23,
  99      >     t46 = t23 * t23, nv = 64)
 100       double precision  xv(nv), t1, t2, t3, t4, an, a1, a2, x1, x2, yy
 101       integer n1, i, j
 102       external randlc
 103       double precision randlc
 104
 105 c---------------------------------------------------------------------
 106 c     Compute the first NV elements of the sequence using RANDLC.
 107 c---------------------------------------------------------------------
 108       t1 = x
 109       n1 = min (n, nv)
 110
 111       do  i = 1, n1
 112          xv(i) = t46 * randlc (t1, a)
 113       enddo
 114
 115 c---------------------------------------------------------------------
 116 c     It is not necessary to compute AN, A1 or A2 unless N is greater than NV.
 117 c---------------------------------------------------------------------
 118       if (n .gt. nv) then
 119
 120 c---------------------------------------------------------------------
 121 c     Compute AN = AA ^ NV (mod 2^46) using successive calls to RANDLC.
 122 c---------------------------------------------------------------------
 123          t1 = a
 124          t2 = r46 * a
 125
 126          do  i = 1, nv - 1
 127             t2 = randlc (t1, a)
 128          enddo
 129
 130          an = t46 * t2
 131
 132 c---------------------------------------------------------------------
 133 c     Break AN into two parts such that AN = 2^23 * A1 + A2.
 134 c---------------------------------------------------------------------
 135          t1 = r23 * an
 136          a1 = aint (t1)
 137          a2 = an - t23 * a1
 138       endif
 139
 140 c---------------------------------------------------------------------
 141 c     Compute N pseudorandom results in batches of size NV.
 142 c---------------------------------------------------------------------
 143       do  j = 0, n - 1, nv
 144          n1 = min (nv, n - j)
 145
 146 c---------------------------------------------------------------------
 147 c     Compute up to NV results based on the current seed vector XV.
 148 c---------------------------------------------------------------------
 149          do  i = 1, n1
 150             y(i+j) = r46 * xv(i)
 151          enddo
 152
 153 c---------------------------------------------------------------------
 154 c     If this is the last pass through the 140 loop, it is not necessary to
 155 c     update the XV vector.
 156 c---------------------------------------------------------------------
 157          if (j + n1 .eq. n) goto 150
 158
 159 c---------------------------------------------------------------------
 160 c     Update the XV vector by multiplying each element by AN (mod 2^46).
 161 c---------------------------------------------------------------------
 162          do  i = 1, nv
 163             t1 = r23 * xv(i)
 164             x1 = aint (t1)
 165             x2 = xv(i) - t23 * x1
 166             t1 = a1 * x2 + a2 * x1
 167             t2 = aint (r23 * t1)
 168             yy = t1 - t23 * t2
 169             t3 = t23 * yy + a2 * x2
 170             t4 = aint (r46 * t3)
 171             xv(i) = t3 - t46 * t4
 172          enddo
 173
 174       enddo
 175
 176 c---------------------------------------------------------------------
 177 c     Save the last seed in X so that subsequent calls to VRANLC will generate
 178 c     a continuous sequence.
 179 c---------------------------------------------------------------------
 180  150  x = xv(n1)
 181
 182       return
 183       end
 184
 185 c----- end of program ------------------------------------------------
 186