ScaLAPACK 2.1  2.1 ScaLAPACK: Scalable Linear Algebra PACKage
dlaran.f
Go to the documentation of this file.
1  DOUBLE PRECISION FUNCTION dlaran( ISEED )
2 *
3 * -- LAPACK auxiliary routine (version 3.1) --
4 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
5 * November 2006
6 *
7 * .. Array Arguments ..
8  INTEGER iseed( 4 )
9 * ..
10 *
11 * Purpose
12 * =======
13 *
14 * DLARAN returns a random real number from a uniform (0,1)
15 * distribution.
16 *
17 * Arguments
18 * =========
19 *
20 * ISEED (input/output) INTEGER array, dimension (4)
21 * On entry, the seed of the random number generator; the array
22 * elements must be between 0 and 4095, and ISEED(4) must be
23 * odd.
24 * On exit, the seed is updated.
25 *
26 * Further Details
27 * ===============
28 *
29 * This routine uses a multiplicative congruential method with modulus
30 * 2**48 and multiplier 33952834046453 (see G.S.Fishman,
31 * 'Multiplicative congruential random number generators with modulus
32 * 2**b: an exhaustive analysis for b = 32 and a partial analysis for
33 * b = 48', Math. Comp. 189, pp 331-344, 1990).
34 *
35 * 48-bit integers are stored in 4 integer array elements with 12 bits
36 * per element. Hence the routine is portable across machines with
37 * integers of 32 bits or more.
38 *
39 * =====================================================================
40 *
41 * .. Parameters ..
42  INTEGER m1, m2, m3, m4
43  parameter( m1 = 494, m2 = 322, m3 = 2508, m4 = 2549 )
44  DOUBLE PRECISION one
45  parameter( one = 1.0d+0 )
46  INTEGER ipw2
47  DOUBLE PRECISION r
48  parameter( ipw2 = 4096, r = one / ipw2 )
49 * ..
50 * .. Local Scalars ..
51  INTEGER it1, it2, it3, it4
52  DOUBLE PRECISION rndout
53 * ..
54 * .. Intrinsic Functions ..
55  INTRINSIC dble, mod
56 * ..
57 * .. Executable Statements ..
58  10 CONTINUE
59 *
60 * multiply the seed by the multiplier modulo 2**48
61 *
62  it4 = iseed( 4 )*m4
63  it3 = it4 / ipw2
64  it4 = it4 - ipw2*it3
65  it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
66  it2 = it3 / ipw2
67  it3 = it3 - ipw2*it2
68  it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
69  it1 = it2 / ipw2
70  it2 = it2 - ipw2*it1
71  it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 +
72  \$ iseed( 4 )*m1
73  it1 = mod( it1, ipw2 )
74 *
75 * return updated seed
76 *
77  iseed( 1 ) = it1
78  iseed( 2 ) = it2
79  iseed( 3 ) = it3
80  iseed( 4 ) = it4
81 *
82 * convert 48-bit integer to a real number in the interval (0,1)
83 *
84  rndout = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
85  \$ ( dble( it4 ) ) ) ) )
86 *
87  IF (rndout.EQ.1.0d+0) THEN
88 * If a real number has n bits of precision, and the first
89 * n bits of the 48-bit integer above happen to be all 1 (which
90 * will occur about once every 2**n calls), then DLARAN will
91 * be rounded to exactly 1.0.
92 * Since DLARAN is not supposed to return exactly 0.0 or 1.0
93 * (and some callers of DLARAN, such as CLARND, depend on that),
94 * the statistically correct thing to do in this situation is
95 * simply to iterate again.
96 * N.B. the case DLARAN = 0.0 should not be possible.
97 *
98  GOTO 10
99  END IF
100 *
101  dlaran = rndout
102  RETURN
103 *
104 * End of DLARAN
105 *
106  END
dlaran
double precision function dlaran(ISEED)
Definition: tools.f:2000