LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
slaruv.f
Go to the documentation of this file.
1 *> \brief \b SLARUV 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 SLARUV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/slaruv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/slaruv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/slaruv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SLARUV( ISEED, N, X )
22 *
23 * .. Scalar Arguments ..
24 * INTEGER N
25 * ..
26 * .. Array Arguments ..
27 * INTEGER ISEED( 4 )
28 * REAL X( N )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> SLARUV 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 SLARNV and CLARNV.
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 REAL 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 *> \date September 2012
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 * =====================================================================
96  SUBROUTINE slaruv( ISEED, N, X )
97 *
98 * -- LAPACK auxiliary routine (version 3.4.2) --
99 * -- LAPACK is a software package provided by Univ. of Tennessee, --
100 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
101 * September 2012
102 *
103 * .. Scalar Arguments ..
104  INTEGER N
105 * ..
106 * .. Array Arguments ..
107  INTEGER ISEED( 4 )
108  REAL X( n )
109 * ..
110 *
111 * =====================================================================
112 *
113 * .. Parameters ..
114  REAL ONE
115  parameter ( one = 1.0e0 )
116  INTEGER LV, IPW2
117  REAL 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 min, mod, real
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*( REAL( it1 )+R*( REAL( it2 )+r*( REAL( it3 )+r*
416  $ REAL( IT4 ) ) ) )
417 *
418  IF (x( i ).EQ.1.0) 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. In IEEE single precision arithmetic,
423 * this will happen relatively often since n = 24.
424 * Since X( I ) is not supposed to return exactly 0.0 or 1.0,
425 * the statistically correct thing to do in this situation is
426 * simply to iterate again.
427 * N.B. the case X( I ) = 0.0 should not be possible.
428  i1 = i1 + 2
429  i2 = i2 + 2
430  i3 = i3 + 2
431  i4 = i4 + 2
432  GOTO 20
433  END IF
434 *
435  10 CONTINUE
436 *
437 * Return final value of seed
438 *
439  iseed( 1 ) = it1
440  iseed( 2 ) = it2
441  iseed( 3 ) = it3
442  iseed( 4 ) = it4
443  RETURN
444 *
445 * End of SLARUV
446 *
447  END
slaruv
subroutine slaruv(ISEED, N, X)
SLARUV returns a vector of n random real numbers from a uniform distribution.
Definition: slaruv.f:97