LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
slaran.f
Go to the documentation of this file.
1 *> \brief \b SLARAN
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * REAL FUNCTION SLARAN( ISEED )
12 *
13 * .. Array Arguments ..
14 * INTEGER ISEED( 4 )
15 * ..
16 *
17 *
18 *> \par Purpose:
19 * =============
20 *>
21 *> \verbatim
22 *>
23 *> SLARAN returns a random real number from a uniform (0,1)
24 *> distribution.
25 *> \endverbatim
26 *
27 * Arguments:
28 * ==========
29 *
30 *> \param[in,out] ISEED
31 *> \verbatim
32 *> ISEED is INTEGER array, dimension (4)
33 *> On entry, the seed of the random number generator; the array
34 *> elements must be between 0 and 4095, and ISEED(4) must be
35 *> odd.
36 *> On exit, the seed is updated.
37 *> \endverbatim
38 *
39 * Authors:
40 * ========
41 *
42 *> \author Univ. of Tennessee
43 *> \author Univ. of California Berkeley
44 *> \author Univ. of Colorado Denver
45 *> \author NAG Ltd.
46 *
47 *> \date November 2011
48 *
49 *> \ingroup real_matgen
50 *
51 *> \par Further Details:
52 * =====================
53 *>
54 *> \verbatim
55 *>
56 *> This routine uses a multiplicative congruential method with modulus
57 *> 2**48 and multiplier 33952834046453 (see G.S.Fishman,
58 *> 'Multiplicative congruential random number generators with modulus
59 *> 2**b: an exhaustive analysis for b = 32 and a partial analysis for
60 *> b = 48', Math. Comp. 189, pp 331-344, 1990).
61 *>
62 *> 48-bit integers are stored in 4 integer array elements with 12 bits
63 *> per element. Hence the routine is portable across machines with
64 *> integers of 32 bits or more.
65 *> \endverbatim
66 *>
67 * =====================================================================
68  REAL FUNCTION slaran( ISEED )
69 *
70 * -- LAPACK auxiliary routine (version 3.4.0) --
71 * -- LAPACK is a software package provided by Univ. of Tennessee, --
72 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
73 * November 2011
74 *
75 * .. Array Arguments ..
76  INTEGER iseed( 4 )
77 * ..
78 *
79 * =====================================================================
80 *
81 * .. Parameters ..
82  INTEGER m1, m2, m3, m4
83  parameter( m1 = 494, m2 = 322, m3 = 2508, m4 = 2549 )
84  REAL one
85  parameter( one = 1.0e+0 )
86  INTEGER ipw2
87  REAL r
88  parameter( ipw2 = 4096, r = one / ipw2 )
89 * ..
90 * .. Local Scalars ..
91  INTEGER it1, it2, it3, it4
92  REAL rndout
93 * ..
94 * .. Intrinsic Functions ..
95  INTRINSIC mod, real
96 * ..
97 * .. Executable Statements ..
98  10 CONTINUE
99 *
100 * multiply the seed by the multiplier modulo 2**48
101 *
102  it4 = iseed( 4 )*m4
103  it3 = it4 / ipw2
104  it4 = it4 - ipw2*it3
105  it3 = it3 + iseed( 3 )*m4 + iseed( 4 )*m3
106  it2 = it3 / ipw2
107  it3 = it3 - ipw2*it2
108  it2 = it2 + iseed( 2 )*m4 + iseed( 3 )*m3 + iseed( 4 )*m2
109  it1 = it2 / ipw2
110  it2 = it2 - ipw2*it1
111  it1 = it1 + iseed( 1 )*m4 + iseed( 2 )*m3 + iseed( 3 )*m2 +
112  $ iseed( 4 )*m1
113  it1 = mod( it1, ipw2 )
114 *
115 * return updated seed
116 *
117  iseed( 1 ) = it1
118  iseed( 2 ) = it2
119  iseed( 3 ) = it3
120  iseed( 4 ) = it4
121 *
122 * convert 48-bit integer to a real number in the interval (0,1)
123 *
124  rndout = r*( REAL( it1 )+r*( REAL( it2 )+r*( REAL( it3 )+r*
125  $ ( REAL( IT4 ) ) ) ) )
126 *
127  IF (rndout.EQ.1.0) THEN
128 * If a real number has n bits of precision, and the first
129 * n bits of the 48-bit integer above happen to be all 1 (which
130 * will occur about once every 2**n calls), then SLARAN will
131 * be rounded to exactly 1.0. In IEEE single precision arithmetic,
132 * this will happen relatively often since n = 24.
133 * Since SLARAN is not supposed to return exactly 0.0 or 1.0
134 * (and some callers of SLARAN, such as CLARND, depend on that),
135 * the statistically correct thing to do in this situation is
136 * simply to iterate again.
137 * N.B. the case SLARAN = 0.0 should not be possible.
138 *
139  goto 10
140  END IF
141 *
142  slaran = rndout
143  RETURN
144 *
145 * End of SLARAN
146 *
147  END