ScaLAPACK 2.1  2.1
ScaLAPACK: Scalable Linear Algebra PACKage
slaran.f
Go to the documentation of this file.
1  REAL FUNCTION SLARAN( 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 * SLARAN 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  REAL one
45  parameter( one = 1.0e+0 )
46  INTEGER ipw2
47  REAL r
48  parameter( ipw2 = 4096, r = one / ipw2 )
49 * ..
50 * .. Local Scalars ..
51  INTEGER it1, it2, it3, it4
52  REAL rndout
53 * ..
54 * .. Intrinsic Functions ..
55  INTRINSIC mod, real
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*( real( it1 )+r*( real( it2 )+r*( real( it3 )+r*
85  $ ( real( it4 ) ) ) ) )
86 *
87  IF (rndout.EQ.1.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 SLARAN will
91 * be rounded to exactly 1.0. In IEEE single precision arithmetic,
92 * this will happen relatively often since n = 24.
93 * Since SLARAN is not supposed to return exactly 0.0 or 1.0
94 * (and some callers of SLARAN, such as CLARND, depend on that),
95 * the statistically correct thing to do in this situation is
96 * simply to iterate again.
97 * N.B. the case SLARAN = 0.0 should not be possible.
98 *
99  GOTO 10
100  END IF
101 *
102  slaran = rndout
103  RETURN
104 *
105 * End of SLARAN
106 *
107  END
slaran
real function slaran(ISEED)
Definition: slaran.f:2