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