LAPACK 3.12.0
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches
dlaruv.f
Go to the documentation of this file.
1*> \brief \b DLARUV returns a vector of n random real numbers from a uniform distribution.
2*
3* =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6* http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DLARUV + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaruv.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaruv.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaruv.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18* Definition:
19* ===========
20*
21* SUBROUTINE DLARUV( ISEED, N, X )
22*
23* .. Scalar Arguments ..
24* INTEGER N
25* ..
26* .. Array Arguments ..
27* INTEGER ISEED( 4 )
28* DOUBLE PRECISION X( N )
29* ..
30*
31*
32*> \par Purpose:
33* =============
34*>
35*> \verbatim
36*>
37*> DLARUV returns a vector of n random real numbers from a uniform (0,1)
38*> distribution (n <= 128).
39*>
40*> This is an auxiliary routine called by DLARNV and ZLARNV.
41*> \endverbatim
42*
43* Arguments:
44* ==========
45*
46*> \param[in,out] ISEED
47*> \verbatim
48*> ISEED is INTEGER array, dimension (4)
49*> On entry, the seed of the random number generator; the array
50*> elements must be between 0 and 4095, and ISEED(4) must be
51*> odd.
52*> On exit, the seed is updated.
53*> \endverbatim
54*>
55*> \param[in] N
56*> \verbatim
57*> N is INTEGER
58*> The number of random numbers to be generated. N <= 128.
59*> \endverbatim
60*>
61*> \param[out] X
62*> \verbatim
63*> X is DOUBLE PRECISION array, dimension (N)
64*> The generated random numbers.
65*> \endverbatim
66*
67* Authors:
68* ========
69*
70*> \author Univ. of Tennessee
71*> \author Univ. of California Berkeley
72*> \author Univ. of Colorado Denver
73*> \author NAG Ltd.
74*
75*> \ingroup laruv
76*
77*> \par Further Details:
78* =====================
79*>
80*> \verbatim
81*>
82*> This routine uses a multiplicative congruential method with modulus
83*> 2**48 and multiplier 33952834046453 (see G.S.Fishman,
84*> 'Multiplicative congruential random number generators with modulus
85*> 2**b: an exhaustive analysis for b = 32 and a partial analysis for
86*> b = 48', Math. Comp. 189, pp 331-344, 1990).
87*>
88*> 48-bit integers are stored in 4 integer array elements with 12 bits
89*> per element. Hence the routine is portable across machines with
90*> integers of 32 bits or more.
91*> \endverbatim
92*>
93* =====================================================================
94 SUBROUTINE dlaruv( ISEED, N, X )
95*
96* -- LAPACK auxiliary routine --
97* -- LAPACK is a software package provided by Univ. of Tennessee, --
98* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
99*
100* .. Scalar Arguments ..
101 INTEGER N
102* ..
103* .. Array Arguments ..
104 INTEGER ISEED( 4 )
105 DOUBLE PRECISION X( N )
106* ..
107*
108* =====================================================================
109*
110* .. Parameters ..
111 DOUBLE PRECISION ONE
112 parameter( one = 1.0d0 )
113 INTEGER LV, IPW2
114 DOUBLE PRECISION R
115 parameter( lv = 128, ipw2 = 4096, r = one / ipw2 )
116* ..
117* .. Local Scalars ..
118 INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J
119* ..
120* .. Local Arrays ..
121 INTEGER MM( LV, 4 )
122* ..
123* .. Intrinsic Functions ..
124 INTRINSIC dble, min, mod
125* ..
126* .. Data statements ..
127 DATA ( mm( 1, j ), j = 1, 4 ) / 494, 322, 2508,
128 $ 2549 /
129 DATA ( mm( 2, j ), j = 1, 4 ) / 2637, 789, 3754,
130 $ 1145 /
131 DATA ( mm( 3, j ), j = 1, 4 ) / 255, 1440, 1766,
132 $ 2253 /
133 DATA ( mm( 4, j ), j = 1, 4 ) / 2008, 752, 3572,
134 $ 305 /
135 DATA ( mm( 5, j ), j = 1, 4 ) / 1253, 2859, 2893,
136 $ 3301 /
137 DATA ( mm( 6, j ), j = 1, 4 ) / 3344, 123, 307,
138 $ 1065 /
139 DATA ( mm( 7, j ), j = 1, 4 ) / 4084, 1848, 1297,
140 $ 3133 /
141 DATA ( mm( 8, j ), j = 1, 4 ) / 1739, 643, 3966,
142 $ 2913 /
143 DATA ( mm( 9, j ), j = 1, 4 ) / 3143, 2405, 758,
144 $ 3285 /
145 DATA ( mm( 10, j ), j = 1, 4 ) / 3468, 2638, 2598,
146 $ 1241 /
147 DATA ( mm( 11, j ), j = 1, 4 ) / 688, 2344, 3406,
148 $ 1197 /
149 DATA ( mm( 12, j ), j = 1, 4 ) / 1657, 46, 2922,
150 $ 3729 /
151 DATA ( mm( 13, j ), j = 1, 4 ) / 1238, 3814, 1038,
152 $ 2501 /
153 DATA ( mm( 14, j ), j = 1, 4 ) / 3166, 913, 2934,
154 $ 1673 /
155 DATA ( mm( 15, j ), j = 1, 4 ) / 1292, 3649, 2091,
156 $ 541 /
157 DATA ( mm( 16, j ), j = 1, 4 ) / 3422, 339, 2451,
158 $ 2753 /
159 DATA ( mm( 17, j ), j = 1, 4 ) / 1270, 3808, 1580,
160 $ 949 /
161 DATA ( mm( 18, j ), j = 1, 4 ) / 2016, 822, 1958,
162 $ 2361 /
163 DATA ( mm( 19, j ), j = 1, 4 ) / 154, 2832, 2055,
164 $ 1165 /
165 DATA ( mm( 20, j ), j = 1, 4 ) / 2862, 3078, 1507,
166 $ 4081 /
167 DATA ( mm( 21, j ), j = 1, 4 ) / 697, 3633, 1078,
168 $ 2725 /
169 DATA ( mm( 22, j ), j = 1, 4 ) / 1706, 2970, 3273,
170 $ 3305 /
171 DATA ( mm( 23, j ), j = 1, 4 ) / 491, 637, 17,
172 $ 3069 /
173 DATA ( mm( 24, j ), j = 1, 4 ) / 931, 2249, 854,
174 $ 3617 /
175 DATA ( mm( 25, j ), j = 1, 4 ) / 1444, 2081, 2916,
176 $ 3733 /
177 DATA ( mm( 26, j ), j = 1, 4 ) / 444, 4019, 3971,
178 $ 409 /
179 DATA ( mm( 27, j ), j = 1, 4 ) / 3577, 1478, 2889,
180 $ 2157 /
181 DATA ( mm( 28, j ), j = 1, 4 ) / 3944, 242, 3831,
182 $ 1361 /
183 DATA ( mm( 29, j ), j = 1, 4 ) / 2184, 481, 2621,
184 $ 3973 /
185 DATA ( mm( 30, j ), j = 1, 4 ) / 1661, 2075, 1541,
186 $ 1865 /
187 DATA ( mm( 31, j ), j = 1, 4 ) / 3482, 4058, 893,
188 $ 2525 /
189 DATA ( mm( 32, j ), j = 1, 4 ) / 657, 622, 736,
190 $ 1409 /
191 DATA ( mm( 33, j ), j = 1, 4 ) / 3023, 3376, 3992,
192 $ 3445 /
193 DATA ( mm( 34, j ), j = 1, 4 ) / 3618, 812, 787,
194 $ 3577 /
195 DATA ( mm( 35, j ), j = 1, 4 ) / 1267, 234, 2125,
196 $ 77 /
197 DATA ( mm( 36, j ), j = 1, 4 ) / 1828, 641, 2364,
198 $ 3761 /
199 DATA ( mm( 37, j ), j = 1, 4 ) / 164, 4005, 2460,
200 $ 2149 /
201 DATA ( mm( 38, j ), j = 1, 4 ) / 3798, 1122, 257,
202 $ 1449 /
203 DATA ( mm( 39, j ), j = 1, 4 ) / 3087, 3135, 1574,
204 $ 3005 /
205 DATA ( mm( 40, j ), j = 1, 4 ) / 2400, 2640, 3912,
206 $ 225 /
207 DATA ( mm( 41, j ), j = 1, 4 ) / 2870, 2302, 1216,
208 $ 85 /
209 DATA ( mm( 42, j ), j = 1, 4 ) / 3876, 40, 3248,
210 $ 3673 /
211 DATA ( mm( 43, j ), j = 1, 4 ) / 1905, 1832, 3401,
212 $ 3117 /
213 DATA ( mm( 44, j ), j = 1, 4 ) / 1593, 2247, 2124,
214 $ 3089 /
215 DATA ( mm( 45, j ), j = 1, 4 ) / 1797, 2034, 2762,
216 $ 1349 /
217 DATA ( mm( 46, j ), j = 1, 4 ) / 1234, 2637, 149,
218 $ 2057 /
219 DATA ( mm( 47, j ), j = 1, 4 ) / 3460, 1287, 2245,
220 $ 413 /
221 DATA ( mm( 48, j ), j = 1, 4 ) / 328, 1691, 166,
222 $ 65 /
223 DATA ( mm( 49, j ), j = 1, 4 ) / 2861, 496, 466,
224 $ 1845 /
225 DATA ( mm( 50, j ), j = 1, 4 ) / 1950, 1597, 4018,
226 $ 697 /
227 DATA ( mm( 51, j ), j = 1, 4 ) / 617, 2394, 1399,
228 $ 3085 /
229 DATA ( mm( 52, j ), j = 1, 4 ) / 2070, 2584, 190,
230 $ 3441 /
231 DATA ( mm( 53, j ), j = 1, 4 ) / 3331, 1843, 2879,
232 $ 1573 /
233 DATA ( mm( 54, j ), j = 1, 4 ) / 769, 336, 153,
234 $ 3689 /
235 DATA ( mm( 55, j ), j = 1, 4 ) / 1558, 1472, 2320,
236 $ 2941 /
237 DATA ( mm( 56, j ), j = 1, 4 ) / 2412, 2407, 18,
238 $ 929 /
239 DATA ( mm( 57, j ), j = 1, 4 ) / 2800, 433, 712,
240 $ 533 /
241 DATA ( mm( 58, j ), j = 1, 4 ) / 189, 2096, 2159,
242 $ 2841 /
243 DATA ( mm( 59, j ), j = 1, 4 ) / 287, 1761, 2318,
244 $ 4077 /
245 DATA ( mm( 60, j ), j = 1, 4 ) / 2045, 2810, 2091,
246 $ 721 /
247 DATA ( mm( 61, j ), j = 1, 4 ) / 1227, 566, 3443,
248 $ 2821 /
249 DATA ( mm( 62, j ), j = 1, 4 ) / 2838, 442, 1510,
250 $ 2249 /
251 DATA ( mm( 63, j ), j = 1, 4 ) / 209, 41, 449,
252 $ 2397 /
253 DATA ( mm( 64, j ), j = 1, 4 ) / 2770, 1238, 1956,
254 $ 2817 /
255 DATA ( mm( 65, j ), j = 1, 4 ) / 3654, 1086, 2201,
256 $ 245 /
257 DATA ( mm( 66, j ), j = 1, 4 ) / 3993, 603, 3137,
258 $ 1913 /
259 DATA ( mm( 67, j ), j = 1, 4 ) / 192, 840, 3399,
260 $ 1997 /
261 DATA ( mm( 68, j ), j = 1, 4 ) / 2253, 3168, 1321,
262 $ 3121 /
263 DATA ( mm( 69, j ), j = 1, 4 ) / 3491, 1499, 2271,
264 $ 997 /
265 DATA ( mm( 70, j ), j = 1, 4 ) / 2889, 1084, 3667,
266 $ 1833 /
267 DATA ( mm( 71, j ), j = 1, 4 ) / 2857, 3438, 2703,
268 $ 2877 /
269 DATA ( mm( 72, j ), j = 1, 4 ) / 2094, 2408, 629,
270 $ 1633 /
271 DATA ( mm( 73, j ), j = 1, 4 ) / 1818, 1589, 2365,
272 $ 981 /
273 DATA ( mm( 74, j ), j = 1, 4 ) / 688, 2391, 2431,
274 $ 2009 /
275 DATA ( mm( 75, j ), j = 1, 4 ) / 1407, 288, 1113,
276 $ 941 /
277 DATA ( mm( 76, j ), j = 1, 4 ) / 634, 26, 3922,
278 $ 2449 /
279 DATA ( mm( 77, j ), j = 1, 4 ) / 3231, 512, 2554,
280 $ 197 /
281 DATA ( mm( 78, j ), j = 1, 4 ) / 815, 1456, 184,
282 $ 2441 /
283 DATA ( mm( 79, j ), j = 1, 4 ) / 3524, 171, 2099,
284 $ 285 /
285 DATA ( mm( 80, j ), j = 1, 4 ) / 1914, 1677, 3228,
286 $ 1473 /
287 DATA ( mm( 81, j ), j = 1, 4 ) / 516, 2657, 4012,
288 $ 2741 /
289 DATA ( mm( 82, j ), j = 1, 4 ) / 164, 2270, 1921,
290 $ 3129 /
291 DATA ( mm( 83, j ), j = 1, 4 ) / 303, 2587, 3452,
292 $ 909 /
293 DATA ( mm( 84, j ), j = 1, 4 ) / 2144, 2961, 3901,
294 $ 2801 /
295 DATA ( mm( 85, j ), j = 1, 4 ) / 3480, 1970, 572,
296 $ 421 /
297 DATA ( mm( 86, j ), j = 1, 4 ) / 119, 1817, 3309,
298 $ 4073 /
299 DATA ( mm( 87, j ), j = 1, 4 ) / 3357, 676, 3171,
300 $ 2813 /
301 DATA ( mm( 88, j ), j = 1, 4 ) / 837, 1410, 817,
302 $ 2337 /
303 DATA ( mm( 89, j ), j = 1, 4 ) / 2826, 3723, 3039,
304 $ 1429 /
305 DATA ( mm( 90, j ), j = 1, 4 ) / 2332, 2803, 1696,
306 $ 1177 /
307 DATA ( mm( 91, j ), j = 1, 4 ) / 2089, 3185, 1256,
308 $ 1901 /
309 DATA ( mm( 92, j ), j = 1, 4 ) / 3780, 184, 3715,
310 $ 81 /
311 DATA ( mm( 93, j ), j = 1, 4 ) / 1700, 663, 2077,
312 $ 1669 /
313 DATA ( mm( 94, j ), j = 1, 4 ) / 3712, 499, 3019,
314 $ 2633 /
315 DATA ( mm( 95, j ), j = 1, 4 ) / 150, 3784, 1497,
316 $ 2269 /
317 DATA ( mm( 96, j ), j = 1, 4 ) / 2000, 1631, 1101,
318 $ 129 /
319 DATA ( mm( 97, j ), j = 1, 4 ) / 3375, 1925, 717,
320 $ 1141 /
321 DATA ( mm( 98, j ), j = 1, 4 ) / 1621, 3912, 51,
322 $ 249 /
323 DATA ( mm( 99, j ), j = 1, 4 ) / 3090, 1398, 981,
324 $ 3917 /
325 DATA ( mm( 100, j ), j = 1, 4 ) / 3765, 1349, 1978,
326 $ 2481 /
327 DATA ( mm( 101, j ), j = 1, 4 ) / 1149, 1441, 1813,
328 $ 3941 /
329 DATA ( mm( 102, j ), j = 1, 4 ) / 3146, 2224, 3881,
330 $ 2217 /
331 DATA ( mm( 103, j ), j = 1, 4 ) / 33, 2411, 76,
332 $ 2749 /
333 DATA ( mm( 104, j ), j = 1, 4 ) / 3082, 1907, 3846,
334 $ 3041 /
335 DATA ( mm( 105, j ), j = 1, 4 ) / 2741, 3192, 3694,
336 $ 1877 /
337 DATA ( mm( 106, j ), j = 1, 4 ) / 359, 2786, 1682,
338 $ 345 /
339 DATA ( mm( 107, j ), j = 1, 4 ) / 3316, 382, 124,
340 $ 2861 /
341 DATA ( mm( 108, j ), j = 1, 4 ) / 1749, 37, 1660,
342 $ 1809 /
343 DATA ( mm( 109, j ), j = 1, 4 ) / 185, 759, 3997,
344 $ 3141 /
345 DATA ( mm( 110, j ), j = 1, 4 ) / 2784, 2948, 479,
346 $ 2825 /
347 DATA ( mm( 111, j ), j = 1, 4 ) / 2202, 1862, 1141,
348 $ 157 /
349 DATA ( mm( 112, j ), j = 1, 4 ) / 2199, 3802, 886,
350 $ 2881 /
351 DATA ( mm( 113, j ), j = 1, 4 ) / 1364, 2423, 3514,
352 $ 3637 /
353 DATA ( mm( 114, j ), j = 1, 4 ) / 1244, 2051, 1301,
354 $ 1465 /
355 DATA ( mm( 115, j ), j = 1, 4 ) / 2020, 2295, 3604,
356 $ 2829 /
357 DATA ( mm( 116, j ), j = 1, 4 ) / 3160, 1332, 1888,
358 $ 2161 /
359 DATA ( mm( 117, j ), j = 1, 4 ) / 2785, 1832, 1836,
360 $ 3365 /
361 DATA ( mm( 118, j ), j = 1, 4 ) / 2772, 2405, 1990,
362 $ 361 /
363 DATA ( mm( 119, j ), j = 1, 4 ) / 1217, 3638, 2058,
364 $ 2685 /
365 DATA ( mm( 120, j ), j = 1, 4 ) / 1822, 3661, 692,
366 $ 3745 /
367 DATA ( mm( 121, j ), j = 1, 4 ) / 1245, 327, 1194,
368 $ 2325 /
369 DATA ( mm( 122, j ), j = 1, 4 ) / 2252, 3660, 20,
370 $ 3609 /
371 DATA ( mm( 123, j ), j = 1, 4 ) / 3904, 716, 3285,
372 $ 3821 /
373 DATA ( mm( 124, j ), j = 1, 4 ) / 2774, 1842, 2046,
374 $ 3537 /
375 DATA ( mm( 125, j ), j = 1, 4 ) / 997, 3987, 2107,
376 $ 517 /
377 DATA ( mm( 126, j ), j = 1, 4 ) / 2573, 1368, 3508,
378 $ 3017 /
379 DATA ( mm( 127, j ), j = 1, 4 ) / 1148, 1848, 3525,
380 $ 2141 /
381 DATA ( mm( 128, j ), j = 1, 4 ) / 545, 2366, 3801,
382 $ 1537 /
383* ..
384* .. Executable Statements ..
385*
386* Quick return for N < 1
387 IF ( n < 1 ) THEN
388 RETURN
389 END IF
390*
391 i1 = iseed( 1 )
392 i2 = iseed( 2 )
393 i3 = iseed( 3 )
394 i4 = iseed( 4 )
395*
396 DO 10 i = 1, min( n, lv )
397*
398 20 CONTINUE
399*
400* Multiply the seed by i-th power of the multiplier modulo 2**48
401*
402 it4 = i4*mm( i, 4 )
403 it3 = it4 / ipw2
404 it4 = it4 - ipw2*it3
405 it3 = it3 + i3*mm( i, 4 ) + i4*mm( i, 3 )
406 it2 = it3 / ipw2
407 it3 = it3 - ipw2*it2
408 it2 = it2 + i2*mm( i, 4 ) + i3*mm( i, 3 ) + i4*mm( i, 2 )
409 it1 = it2 / ipw2
410 it2 = it2 - ipw2*it1
411 it1 = it1 + i1*mm( i, 4 ) + i2*mm( i, 3 ) + i3*mm( i, 2 ) +
412 $ i4*mm( i, 1 )
413 it1 = mod( it1, ipw2 )
414*
415* Convert 48-bit integer to a real number in the interval (0,1)
416*
417 x( i ) = r*( dble( it1 )+r*( dble( it2 )+r*( dble( it3 )+r*
418 $ dble( it4 ) ) ) )
419*
420 IF (x( i ).EQ.1.0d0) THEN
421* If a real number has n bits of precision, and the first
422* n bits of the 48-bit integer above happen to be all 1 (which
423* will occur about once every 2**n calls), then X( I ) will
424* be rounded to exactly 1.0.
425* Since X( I ) is not supposed to return exactly 0.0 or 1.0,
426* the statistically correct thing to do in this situation is
427* simply to iterate again.
428* N.B. the case X( I ) = 0.0 should not be possible.
429 i1 = i1 + 2
430 i2 = i2 + 2
431 i3 = i3 + 2
432 i4 = i4 + 2
433 GOTO 20
434 END IF
435*
436 10 CONTINUE
437*
438* Return final value of seed
439*
440 iseed( 1 ) = it1
441 iseed( 2 ) = it2
442 iseed( 3 ) = it3
443 iseed( 4 ) = it4
444 RETURN
445*
446* End of DLARUV
447*
448 END
dlaruv
subroutine dlaruv(iseed, n, x)
DLARUV returns a vector of n random real numbers from a uniform distribution.
Definition dlaruv.f:95