examples/smpi/NAS/common/randdp.f

   1 c---------------------------------------------------------------------
   2 c---------------------------------------------------------------------
   3
   4       double precision function randlc (x, a)
   5
   6 c---------------------------------------------------------------------
   7 c---------------------------------------------------------------------
   8
   9 c---------------------------------------------------------------------
  10 c
  11 c   This routine returns a uniform pseudorandom double precision number in the
  12 c   range (0, 1) by using the linear congruential generator
  13 c
  14 c   x_{k+1} = a x_k  (mod 2^46)
  15 c
  16 c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  17 c   before repeating.  The argument A is the same as 'a' in the above formula,
  18 c   and X is the same as x_0.  A and X must be odd double precision integers
  19 c   in the range (1, 2^46).  The returned value RANDLC is normalized to be
  20 c   between 0 and 1, i.e. RANDLC = 2^(-46) * x_1.  X is updated to contain
  21 c   the new seed x_1, so that subsequent calls to RANDLC using the same
  22 c   arguments will generate a continuous sequence.
  23 c
  24 c   This routine should produce the same results on any computer with at least
  25 c   48 mantissa bits in double precision floating point data.  On 64 bit
  26 c   systems, double precision should be disabled.
  27 c
  28 c   David H. Bailey     October 26, 1990
  29 c
  30 c---------------------------------------------------------------------
  31
  32       implicit none
  33
  34       double precision r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
  35       parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
  36      >  t46 = t23 ** 2)
  37
  38 c---------------------------------------------------------------------
  39 c   Break A into two parts such that A = 2^23 * A1 + A2.
  40 c---------------------------------------------------------------------
  41       t1 = r23 * a
  42       a1 = int (t1)
  43       a2 = a - t23 * a1
  44
  45 c---------------------------------------------------------------------
  46 c   Break X into two parts such that X = 2^23 * X1 + X2, compute
  47 c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
  48 c   X = 2^23 * Z + A2 * X2  (mod 2^46).
  49 c---------------------------------------------------------------------
  50       t1 = r23 * x
  51       x1 = int (t1)
  52       x2 = x - t23 * x1
  53       t1 = a1 * x2 + a2 * x1
  54       t2 = int (r23 * t1)
  55       z = t1 - t23 * t2
  56       t3 = t23 * z + a2 * x2
  57       t4 = int (r46 * t3)
  58       x = t3 - t46 * t4
  59       randlc = r46 * x
  60
  61       return
  62       end
  63
  64
  65
  66
  67 c---------------------------------------------------------------------
  68 c---------------------------------------------------------------------
  69
  70       subroutine vranlc (n, x, a, y)
  71
  72 c---------------------------------------------------------------------
  73 c---------------------------------------------------------------------
  74
  75 c---------------------------------------------------------------------
  76 c
  77 c   This routine generates N uniform pseudorandom double precision numbers in
  78 c   the range (0, 1) by using the linear congruential generator
  79 c
  80 c   x_{k+1} = a x_k  (mod 2^46)
  81 c
  82 c   where 0 < x_k < 2^46 and 0 < a < 2^46.  This scheme generates 2^44 numbers
  83 c   before repeating.  The argument A is the same as 'a' in the above formula,
  84 c   and X is the same as x_0.  A and X must be odd double precision integers
  85 c   in the range (1, 2^46).  The N results are placed in Y and are normalized
  86 c   to be between 0 and 1.  X is updated to contain the new seed, so that
  87 c   subsequent calls to VRANLC using the same arguments will generate a
  88 c   continuous sequence.  If N is zero, only initialization is performed, and
  89 c   the variables X, A and Y are ignored.
  90 c
  91 c   This routine is the standard version designed for scalar or RISC systems.
  92 c   However, it should produce the same results on any single processor
  93 c   computer with at least 48 mantissa bits in double precision floating point
  94 c   data.  On 64 bit systems, double precision should be disabled.
  95 c
  96 c---------------------------------------------------------------------
  97
  98       implicit none
  99
 100       integer i,n
 101       double precision y,r23,r46,t23,t46,a,x,t1,t2,t3,t4,a1,a2,x1,x2,z
 102       dimension y(*)
 103       parameter (r23 = 0.5d0 ** 23, r46 = r23 ** 2, t23 = 2.d0 ** 23,
 104      >  t46 = t23 ** 2)
 105
 106
 107 c---------------------------------------------------------------------
 108 c   Break A into two parts such that A = 2^23 * A1 + A2.
 109 c---------------------------------------------------------------------
 110       t1 = r23 * a
 111       a1 = int (t1)
 112       a2 = a - t23 * a1
 113
 114 c---------------------------------------------------------------------
 115 c   Generate N results.   This loop is not vectorizable.
 116 c---------------------------------------------------------------------
 117       do i = 1, n
 118
 119 c---------------------------------------------------------------------
 120 c   Break X into two parts such that X = 2^23 * X1 + X2, compute
 121 c   Z = A1 * X2 + A2 * X1  (mod 2^23), and then
 122 c   X = 2^23 * Z + A2 * X2  (mod 2^46).
 123 c---------------------------------------------------------------------
 124         t1 = r23 * x
 125         x1 = int (t1)
 126         x2 = x - t23 * x1
 127         t1 = a1 * x2 + a2 * x1
 128         t2 = int (r23 * t1)
 129         z = t1 - t23 * t2
 130         t3 = t23 * z + a2 * x2
 131         t4 = int (r46 * t3)
 132         x = t3 - t46 * t4
 133         y(i) = r46 * x
 134       enddo
 135
 136       return
 137       end