![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLARUV 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DLARUV + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaruv.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaruv.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaruv.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DLARUV( ISEED, N, X ) 00022 * 00023 * .. Scalar Arguments .. 00024 * INTEGER N 00025 * .. 00026 * .. Array Arguments .. 00027 * INTEGER ISEED( 4 ) 00028 * DOUBLE PRECISION X( N ) 00029 * .. 00030 * 00031 * 00032 *> \par Purpose: 00033 * ============= 00034 *> 00035 *> \verbatim 00036 *> 00037 *> DLARUV returns a vector of n random real numbers from a uniform (0,1) 00038 *> distribution (n <= 128). 00039 *> 00040 *> This is an auxiliary routine called by DLARNV and ZLARNV. 00041 *> \endverbatim 00042 * 00043 * Arguments: 00044 * ========== 00045 * 00046 *> \param[in,out] ISEED 00047 *> \verbatim 00048 *> ISEED is INTEGER array, dimension (4) 00049 *> On entry, the seed of the random number generator; the array 00050 *> elements must be between 0 and 4095, and ISEED(4) must be 00051 *> odd. 00052 *> On exit, the seed is updated. 00053 *> \endverbatim 00054 *> 00055 *> \param[in] N 00056 *> \verbatim 00057 *> N is INTEGER 00058 *> The number of random numbers to be generated. N <= 128. 00059 *> \endverbatim 00060 *> 00061 *> \param[out] X 00062 *> \verbatim 00063 *> X is DOUBLE PRECISION array, dimension (N) 00064 *> The generated random numbers. 00065 *> \endverbatim 00066 * 00067 * Authors: 00068 * ======== 00069 * 00070 *> \author Univ. of Tennessee 00071 *> \author Univ. of California Berkeley 00072 *> \author Univ. of Colorado Denver 00073 *> \author NAG Ltd. 00074 * 00075 *> \date November 2011 00076 * 00077 *> \ingroup auxOTHERauxiliary 00078 * 00079 *> \par Further Details: 00080 * ===================== 00081 *> 00082 *> \verbatim 00083 *> 00084 *> This routine uses a multiplicative congruential method with modulus 00085 *> 2**48 and multiplier 33952834046453 (see G.S.Fishman, 00086 *> 'Multiplicative congruential random number generators with modulus 00087 *> 2**b: an exhaustive analysis for b = 32 and a partial analysis for 00088 *> b = 48', Math. Comp. 189, pp 331-344, 1990). 00089 *> 00090 *> 48-bit integers are stored in 4 integer array elements with 12 bits 00091 *> per element. Hence the routine is portable across machines with 00092 *> integers of 32 bits or more. 00093 *> \endverbatim 00094 *> 00095 * ===================================================================== 00096 SUBROUTINE DLARUV( ISEED, N, X ) 00097 * 00098 * -- LAPACK auxiliary routine (version 3.4.0) -- 00099 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00100 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00101 * November 2011 00102 * 00103 * .. Scalar Arguments .. 00104 INTEGER N 00105 * .. 00106 * .. Array Arguments .. 00107 INTEGER ISEED( 4 ) 00108 DOUBLE PRECISION X( N ) 00109 * .. 00110 * 00111 * ===================================================================== 00112 * 00113 * .. Parameters .. 00114 DOUBLE PRECISION ONE 00115 PARAMETER ( ONE = 1.0D0 ) 00116 INTEGER LV, IPW2 00117 DOUBLE PRECISION R 00118 PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 ) 00119 * .. 00120 * .. Local Scalars .. 00121 INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J 00122 * .. 00123 * .. Local Arrays .. 00124 INTEGER MM( LV, 4 ) 00125 * .. 00126 * .. Intrinsic Functions .. 00127 INTRINSIC DBLE, MIN, MOD 00128 * .. 00129 * .. Data statements .. 00130 DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508, 00131 $ 2549 / 00132 DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754, 00133 $ 1145 / 00134 DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766, 00135 $ 2253 / 00136 DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572, 00137 $ 305 / 00138 DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893, 00139 $ 3301 / 00140 DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307, 00141 $ 1065 / 00142 DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297, 00143 $ 3133 / 00144 DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966, 00145 $ 2913 / 00146 DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758, 00147 $ 3285 / 00148 DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598, 00149 $ 1241 / 00150 DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406, 00151 $ 1197 / 00152 DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922, 00153 $ 3729 / 00154 DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038, 00155 $ 2501 / 00156 DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934, 00157 $ 1673 / 00158 DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091, 00159 $ 541 / 00160 DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451, 00161 $ 2753 / 00162 DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580, 00163 $ 949 / 00164 DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958, 00165 $ 2361 / 00166 DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055, 00167 $ 1165 / 00168 DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507, 00169 $ 4081 / 00170 DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078, 00171 $ 2725 / 00172 DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273, 00173 $ 3305 / 00174 DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17, 00175 $ 3069 / 00176 DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854, 00177 $ 3617 / 00178 DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916, 00179 $ 3733 / 00180 DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971, 00181 $ 409 / 00182 DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889, 00183 $ 2157 / 00184 DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831, 00185 $ 1361 / 00186 DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621, 00187 $ 3973 / 00188 DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541, 00189 $ 1865 / 00190 DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893, 00191 $ 2525 / 00192 DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736, 00193 $ 1409 / 00194 DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992, 00195 $ 3445 / 00196 DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787, 00197 $ 3577 / 00198 DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125, 00199 $ 77 / 00200 DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364, 00201 $ 3761 / 00202 DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460, 00203 $ 2149 / 00204 DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257, 00205 $ 1449 / 00206 DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574, 00207 $ 3005 / 00208 DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912, 00209 $ 225 / 00210 DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216, 00211 $ 85 / 00212 DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248, 00213 $ 3673 / 00214 DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401, 00215 $ 3117 / 00216 DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124, 00217 $ 3089 / 00218 DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762, 00219 $ 1349 / 00220 DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149, 00221 $ 2057 / 00222 DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245, 00223 $ 413 / 00224 DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166, 00225 $ 65 / 00226 DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466, 00227 $ 1845 / 00228 DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018, 00229 $ 697 / 00230 DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399, 00231 $ 3085 / 00232 DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190, 00233 $ 3441 / 00234 DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879, 00235 $ 1573 / 00236 DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153, 00237 $ 3689 / 00238 DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320, 00239 $ 2941 / 00240 DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18, 00241 $ 929 / 00242 DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712, 00243 $ 533 / 00244 DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159, 00245 $ 2841 / 00246 DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318, 00247 $ 4077 / 00248 DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091, 00249 $ 721 / 00250 DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443, 00251 $ 2821 / 00252 DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510, 00253 $ 2249 / 00254 DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449, 00255 $ 2397 / 00256 DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956, 00257 $ 2817 / 00258 DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201, 00259 $ 245 / 00260 DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137, 00261 $ 1913 / 00262 DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399, 00263 $ 1997 / 00264 DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321, 00265 $ 3121 / 00266 DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271, 00267 $ 997 / 00268 DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667, 00269 $ 1833 / 00270 DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703, 00271 $ 2877 / 00272 DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629, 00273 $ 1633 / 00274 DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365, 00275 $ 981 / 00276 DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431, 00277 $ 2009 / 00278 DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113, 00279 $ 941 / 00280 DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922, 00281 $ 2449 / 00282 DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554, 00283 $ 197 / 00284 DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184, 00285 $ 2441 / 00286 DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099, 00287 $ 285 / 00288 DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228, 00289 $ 1473 / 00290 DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012, 00291 $ 2741 / 00292 DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921, 00293 $ 3129 / 00294 DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452, 00295 $ 909 / 00296 DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901, 00297 $ 2801 / 00298 DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572, 00299 $ 421 / 00300 DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309, 00301 $ 4073 / 00302 DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171, 00303 $ 2813 / 00304 DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817, 00305 $ 2337 / 00306 DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039, 00307 $ 1429 / 00308 DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696, 00309 $ 1177 / 00310 DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256, 00311 $ 1901 / 00312 DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715, 00313 $ 81 / 00314 DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077, 00315 $ 1669 / 00316 DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019, 00317 $ 2633 / 00318 DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497, 00319 $ 2269 / 00320 DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101, 00321 $ 129 / 00322 DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717, 00323 $ 1141 / 00324 DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51, 00325 $ 249 / 00326 DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981, 00327 $ 3917 / 00328 DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978, 00329 $ 2481 / 00330 DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813, 00331 $ 3941 / 00332 DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881, 00333 $ 2217 / 00334 DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76, 00335 $ 2749 / 00336 DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846, 00337 $ 3041 / 00338 DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694, 00339 $ 1877 / 00340 DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682, 00341 $ 345 / 00342 DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124, 00343 $ 2861 / 00344 DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660, 00345 $ 1809 / 00346 DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997, 00347 $ 3141 / 00348 DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479, 00349 $ 2825 / 00350 DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141, 00351 $ 157 / 00352 DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886, 00353 $ 2881 / 00354 DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514, 00355 $ 3637 / 00356 DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301, 00357 $ 1465 / 00358 DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604, 00359 $ 2829 / 00360 DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888, 00361 $ 2161 / 00362 DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836, 00363 $ 3365 / 00364 DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990, 00365 $ 361 / 00366 DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058, 00367 $ 2685 / 00368 DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692, 00369 $ 3745 / 00370 DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194, 00371 $ 2325 / 00372 DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20, 00373 $ 3609 / 00374 DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285, 00375 $ 3821 / 00376 DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046, 00377 $ 3537 / 00378 DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107, 00379 $ 517 / 00380 DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508, 00381 $ 3017 / 00382 DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525, 00383 $ 2141 / 00384 DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801, 00385 $ 1537 / 00386 * .. 00387 * .. Executable Statements .. 00388 * 00389 I1 = ISEED( 1 ) 00390 I2 = ISEED( 2 ) 00391 I3 = ISEED( 3 ) 00392 I4 = ISEED( 4 ) 00393 * 00394 DO 10 I = 1, MIN( N, LV ) 00395 * 00396 20 CONTINUE 00397 * 00398 * Multiply the seed by i-th power of the multiplier modulo 2**48 00399 * 00400 IT4 = I4*MM( I, 4 ) 00401 IT3 = IT4 / IPW2 00402 IT4 = IT4 - IPW2*IT3 00403 IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 ) 00404 IT2 = IT3 / IPW2 00405 IT3 = IT3 - IPW2*IT2 00406 IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 ) 00407 IT1 = IT2 / IPW2 00408 IT2 = IT2 - IPW2*IT1 00409 IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) + 00410 $ I4*MM( I, 1 ) 00411 IT1 = MOD( IT1, IPW2 ) 00412 * 00413 * Convert 48-bit integer to a real number in the interval (0,1) 00414 * 00415 X( I ) = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* 00416 $ DBLE( IT4 ) ) ) ) 00417 * 00418 IF (X( I ).EQ.1.0D0) THEN 00419 * If a real number has n bits of precision, and the first 00420 * n bits of the 48-bit integer above happen to be all 1 (which 00421 * will occur about once every 2**n calls), then X( I ) will 00422 * be rounded to exactly 1.0. 00423 * Since X( I ) is not supposed to return exactly 0.0 or 1.0, 00424 * the statistically correct thing to do in this situation is 00425 * simply to iterate again. 00426 * N.B. the case X( I ) = 0.0 should not be possible. 00427 I1 = I1 + 2 00428 I2 = I2 + 2 00429 I3 = I3 + 2 00430 I4 = I4 + 2 00431 GOTO 20 00432 END IF 00433 * 00434 10 CONTINUE 00435 * 00436 * Return final value of seed 00437 * 00438 ISEED( 1 ) = IT1 00439 ISEED( 2 ) = IT2 00440 ISEED( 3 ) = IT3 00441 ISEED( 4 ) = IT4 00442 RETURN 00443 * 00444 * End of DLARUV 00445 * 00446 END