SUBROUTINE DLARUV( ISEED, N, X ) * * -- LAPACK auxiliary routine (version 3.0) -- * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., * Courant Institute, Argonne National Lab, and Rice University * October 31, 1992 * * .. Scalar Arguments .. INTEGER N * .. * .. Array Arguments .. INTEGER ISEED( 4 ) DOUBLE PRECISION X( N ) * .. * * Purpose * ======= * * DLARUV returns a vector of n random real numbers from a uniform (0,1) * distribution (n <= 128). * * This is an auxiliary routine called by DLARNV and ZLARNV. * * Arguments * ========= * * ISEED (input/output) INTEGER array, dimension (4) * On entry, the seed of the random number generator; the array * elements must be between 0 and 4095, and ISEED(4) must be * odd. * On exit, the seed is updated. * * N (input) INTEGER * The number of random numbers to be generated. N <= 128. * * X (output) DOUBLE PRECISION array, dimension (N) * The generated random numbers. * * Further Details * =============== * * This routine uses a multiplicative congruential method with modulus * 2**48 and multiplier 33952834046453 (see G.S.Fishman, * 'Multiplicative congruential random number generators with modulus * 2**b: an exhaustive analysis for b = 32 and a partial analysis for * b = 48', Math. Comp. 189, pp 331-344, 1990). * * 48-bit integers are stored in 4 integer array elements with 12 bits * per element. Hence the routine is portable across machines with * integers of 32 bits or more. * * ===================================================================== * * .. Parameters .. DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D0 ) INTEGER LV, IPW2 DOUBLE PRECISION R PARAMETER ( LV = 128, IPW2 = 4096, R = ONE / IPW2 ) * .. * .. Local Scalars .. INTEGER I, I1, I2, I3, I4, IT1, IT2, IT3, IT4, J * .. * .. Local Arrays .. INTEGER MM( LV, 4 ) * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MIN, MOD * .. * .. Data statements .. DATA ( MM( 1, J ), J = 1, 4 ) / 494, 322, 2508, \$ 2549 / DATA ( MM( 2, J ), J = 1, 4 ) / 2637, 789, 3754, \$ 1145 / DATA ( MM( 3, J ), J = 1, 4 ) / 255, 1440, 1766, \$ 2253 / DATA ( MM( 4, J ), J = 1, 4 ) / 2008, 752, 3572, \$ 305 / DATA ( MM( 5, J ), J = 1, 4 ) / 1253, 2859, 2893, \$ 3301 / DATA ( MM( 6, J ), J = 1, 4 ) / 3344, 123, 307, \$ 1065 / DATA ( MM( 7, J ), J = 1, 4 ) / 4084, 1848, 1297, \$ 3133 / DATA ( MM( 8, J ), J = 1, 4 ) / 1739, 643, 3966, \$ 2913 / DATA ( MM( 9, J ), J = 1, 4 ) / 3143, 2405, 758, \$ 3285 / DATA ( MM( 10, J ), J = 1, 4 ) / 3468, 2638, 2598, \$ 1241 / DATA ( MM( 11, J ), J = 1, 4 ) / 688, 2344, 3406, \$ 1197 / DATA ( MM( 12, J ), J = 1, 4 ) / 1657, 46, 2922, \$ 3729 / DATA ( MM( 13, J ), J = 1, 4 ) / 1238, 3814, 1038, \$ 2501 / DATA ( MM( 14, J ), J = 1, 4 ) / 3166, 913, 2934, \$ 1673 / DATA ( MM( 15, J ), J = 1, 4 ) / 1292, 3649, 2091, \$ 541 / DATA ( MM( 16, J ), J = 1, 4 ) / 3422, 339, 2451, \$ 2753 / DATA ( MM( 17, J ), J = 1, 4 ) / 1270, 3808, 1580, \$ 949 / DATA ( MM( 18, J ), J = 1, 4 ) / 2016, 822, 1958, \$ 2361 / DATA ( MM( 19, J ), J = 1, 4 ) / 154, 2832, 2055, \$ 1165 / DATA ( MM( 20, J ), J = 1, 4 ) / 2862, 3078, 1507, \$ 4081 / DATA ( MM( 21, J ), J = 1, 4 ) / 697, 3633, 1078, \$ 2725 / DATA ( MM( 22, J ), J = 1, 4 ) / 1706, 2970, 3273, \$ 3305 / DATA ( MM( 23, J ), J = 1, 4 ) / 491, 637, 17, \$ 3069 / DATA ( MM( 24, J ), J = 1, 4 ) / 931, 2249, 854, \$ 3617 / DATA ( MM( 25, J ), J = 1, 4 ) / 1444, 2081, 2916, \$ 3733 / DATA ( MM( 26, J ), J = 1, 4 ) / 444, 4019, 3971, \$ 409 / DATA ( MM( 27, J ), J = 1, 4 ) / 3577, 1478, 2889, \$ 2157 / DATA ( MM( 28, J ), J = 1, 4 ) / 3944, 242, 3831, \$ 1361 / DATA ( MM( 29, J ), J = 1, 4 ) / 2184, 481, 2621, \$ 3973 / DATA ( MM( 30, J ), J = 1, 4 ) / 1661, 2075, 1541, \$ 1865 / DATA ( MM( 31, J ), J = 1, 4 ) / 3482, 4058, 893, \$ 2525 / DATA ( MM( 32, J ), J = 1, 4 ) / 657, 622, 736, \$ 1409 / DATA ( MM( 33, J ), J = 1, 4 ) / 3023, 3376, 3992, \$ 3445 / DATA ( MM( 34, J ), J = 1, 4 ) / 3618, 812, 787, \$ 3577 / DATA ( MM( 35, J ), J = 1, 4 ) / 1267, 234, 2125, \$ 77 / DATA ( MM( 36, J ), J = 1, 4 ) / 1828, 641, 2364, \$ 3761 / DATA ( MM( 37, J ), J = 1, 4 ) / 164, 4005, 2460, \$ 2149 / DATA ( MM( 38, J ), J = 1, 4 ) / 3798, 1122, 257, \$ 1449 / DATA ( MM( 39, J ), J = 1, 4 ) / 3087, 3135, 1574, \$ 3005 / DATA ( MM( 40, J ), J = 1, 4 ) / 2400, 2640, 3912, \$ 225 / DATA ( MM( 41, J ), J = 1, 4 ) / 2870, 2302, 1216, \$ 85 / DATA ( MM( 42, J ), J = 1, 4 ) / 3876, 40, 3248, \$ 3673 / DATA ( MM( 43, J ), J = 1, 4 ) / 1905, 1832, 3401, \$ 3117 / DATA ( MM( 44, J ), J = 1, 4 ) / 1593, 2247, 2124, \$ 3089 / DATA ( MM( 45, J ), J = 1, 4 ) / 1797, 2034, 2762, \$ 1349 / DATA ( MM( 46, J ), J = 1, 4 ) / 1234, 2637, 149, \$ 2057 / DATA ( MM( 47, J ), J = 1, 4 ) / 3460, 1287, 2245, \$ 413 / DATA ( MM( 48, J ), J = 1, 4 ) / 328, 1691, 166, \$ 65 / DATA ( MM( 49, J ), J = 1, 4 ) / 2861, 496, 466, \$ 1845 / DATA ( MM( 50, J ), J = 1, 4 ) / 1950, 1597, 4018, \$ 697 / DATA ( MM( 51, J ), J = 1, 4 ) / 617, 2394, 1399, \$ 3085 / DATA ( MM( 52, J ), J = 1, 4 ) / 2070, 2584, 190, \$ 3441 / DATA ( MM( 53, J ), J = 1, 4 ) / 3331, 1843, 2879, \$ 1573 / DATA ( MM( 54, J ), J = 1, 4 ) / 769, 336, 153, \$ 3689 / DATA ( MM( 55, J ), J = 1, 4 ) / 1558, 1472, 2320, \$ 2941 / DATA ( MM( 56, J ), J = 1, 4 ) / 2412, 2407, 18, \$ 929 / DATA ( MM( 57, J ), J = 1, 4 ) / 2800, 433, 712, \$ 533 / DATA ( MM( 58, J ), J = 1, 4 ) / 189, 2096, 2159, \$ 2841 / DATA ( MM( 59, J ), J = 1, 4 ) / 287, 1761, 2318, \$ 4077 / DATA ( MM( 60, J ), J = 1, 4 ) / 2045, 2810, 2091, \$ 721 / DATA ( MM( 61, J ), J = 1, 4 ) / 1227, 566, 3443, \$ 2821 / DATA ( MM( 62, J ), J = 1, 4 ) / 2838, 442, 1510, \$ 2249 / DATA ( MM( 63, J ), J = 1, 4 ) / 209, 41, 449, \$ 2397 / DATA ( MM( 64, J ), J = 1, 4 ) / 2770, 1238, 1956, \$ 2817 / DATA ( MM( 65, J ), J = 1, 4 ) / 3654, 1086, 2201, \$ 245 / DATA ( MM( 66, J ), J = 1, 4 ) / 3993, 603, 3137, \$ 1913 / DATA ( MM( 67, J ), J = 1, 4 ) / 192, 840, 3399, \$ 1997 / DATA ( MM( 68, J ), J = 1, 4 ) / 2253, 3168, 1321, \$ 3121 / DATA ( MM( 69, J ), J = 1, 4 ) / 3491, 1499, 2271, \$ 997 / DATA ( MM( 70, J ), J = 1, 4 ) / 2889, 1084, 3667, \$ 1833 / DATA ( MM( 71, J ), J = 1, 4 ) / 2857, 3438, 2703, \$ 2877 / DATA ( MM( 72, J ), J = 1, 4 ) / 2094, 2408, 629, \$ 1633 / DATA ( MM( 73, J ), J = 1, 4 ) / 1818, 1589, 2365, \$ 981 / DATA ( MM( 74, J ), J = 1, 4 ) / 688, 2391, 2431, \$ 2009 / DATA ( MM( 75, J ), J = 1, 4 ) / 1407, 288, 1113, \$ 941 / DATA ( MM( 76, J ), J = 1, 4 ) / 634, 26, 3922, \$ 2449 / DATA ( MM( 77, J ), J = 1, 4 ) / 3231, 512, 2554, \$ 197 / DATA ( MM( 78, J ), J = 1, 4 ) / 815, 1456, 184, \$ 2441 / DATA ( MM( 79, J ), J = 1, 4 ) / 3524, 171, 2099, \$ 285 / DATA ( MM( 80, J ), J = 1, 4 ) / 1914, 1677, 3228, \$ 1473 / DATA ( MM( 81, J ), J = 1, 4 ) / 516, 2657, 4012, \$ 2741 / DATA ( MM( 82, J ), J = 1, 4 ) / 164, 2270, 1921, \$ 3129 / DATA ( MM( 83, J ), J = 1, 4 ) / 303, 2587, 3452, \$ 909 / DATA ( MM( 84, J ), J = 1, 4 ) / 2144, 2961, 3901, \$ 2801 / DATA ( MM( 85, J ), J = 1, 4 ) / 3480, 1970, 572, \$ 421 / DATA ( MM( 86, J ), J = 1, 4 ) / 119, 1817, 3309, \$ 4073 / DATA ( MM( 87, J ), J = 1, 4 ) / 3357, 676, 3171, \$ 2813 / DATA ( MM( 88, J ), J = 1, 4 ) / 837, 1410, 817, \$ 2337 / DATA ( MM( 89, J ), J = 1, 4 ) / 2826, 3723, 3039, \$ 1429 / DATA ( MM( 90, J ), J = 1, 4 ) / 2332, 2803, 1696, \$ 1177 / DATA ( MM( 91, J ), J = 1, 4 ) / 2089, 3185, 1256, \$ 1901 / DATA ( MM( 92, J ), J = 1, 4 ) / 3780, 184, 3715, \$ 81 / DATA ( MM( 93, J ), J = 1, 4 ) / 1700, 663, 2077, \$ 1669 / DATA ( MM( 94, J ), J = 1, 4 ) / 3712, 499, 3019, \$ 2633 / DATA ( MM( 95, J ), J = 1, 4 ) / 150, 3784, 1497, \$ 2269 / DATA ( MM( 96, J ), J = 1, 4 ) / 2000, 1631, 1101, \$ 129 / DATA ( MM( 97, J ), J = 1, 4 ) / 3375, 1925, 717, \$ 1141 / DATA ( MM( 98, J ), J = 1, 4 ) / 1621, 3912, 51, \$ 249 / DATA ( MM( 99, J ), J = 1, 4 ) / 3090, 1398, 981, \$ 3917 / DATA ( MM( 100, J ), J = 1, 4 ) / 3765, 1349, 1978, \$ 2481 / DATA ( MM( 101, J ), J = 1, 4 ) / 1149, 1441, 1813, \$ 3941 / DATA ( MM( 102, J ), J = 1, 4 ) / 3146, 2224, 3881, \$ 2217 / DATA ( MM( 103, J ), J = 1, 4 ) / 33, 2411, 76, \$ 2749 / DATA ( MM( 104, J ), J = 1, 4 ) / 3082, 1907, 3846, \$ 3041 / DATA ( MM( 105, J ), J = 1, 4 ) / 2741, 3192, 3694, \$ 1877 / DATA ( MM( 106, J ), J = 1, 4 ) / 359, 2786, 1682, \$ 345 / DATA ( MM( 107, J ), J = 1, 4 ) / 3316, 382, 124, \$ 2861 / DATA ( MM( 108, J ), J = 1, 4 ) / 1749, 37, 1660, \$ 1809 / DATA ( MM( 109, J ), J = 1, 4 ) / 185, 759, 3997, \$ 3141 / DATA ( MM( 110, J ), J = 1, 4 ) / 2784, 2948, 479, \$ 2825 / DATA ( MM( 111, J ), J = 1, 4 ) / 2202, 1862, 1141, \$ 157 / DATA ( MM( 112, J ), J = 1, 4 ) / 2199, 3802, 886, \$ 2881 / DATA ( MM( 113, J ), J = 1, 4 ) / 1364, 2423, 3514, \$ 3637 / DATA ( MM( 114, J ), J = 1, 4 ) / 1244, 2051, 1301, \$ 1465 / DATA ( MM( 115, J ), J = 1, 4 ) / 2020, 2295, 3604, \$ 2829 / DATA ( MM( 116, J ), J = 1, 4 ) / 3160, 1332, 1888, \$ 2161 / DATA ( MM( 117, J ), J = 1, 4 ) / 2785, 1832, 1836, \$ 3365 / DATA ( MM( 118, J ), J = 1, 4 ) / 2772, 2405, 1990, \$ 361 / DATA ( MM( 119, J ), J = 1, 4 ) / 1217, 3638, 2058, \$ 2685 / DATA ( MM( 120, J ), J = 1, 4 ) / 1822, 3661, 692, \$ 3745 / DATA ( MM( 121, J ), J = 1, 4 ) / 1245, 327, 1194, \$ 2325 / DATA ( MM( 122, J ), J = 1, 4 ) / 2252, 3660, 20, \$ 3609 / DATA ( MM( 123, J ), J = 1, 4 ) / 3904, 716, 3285, \$ 3821 / DATA ( MM( 124, J ), J = 1, 4 ) / 2774, 1842, 2046, \$ 3537 / DATA ( MM( 125, J ), J = 1, 4 ) / 997, 3987, 2107, \$ 517 / DATA ( MM( 126, J ), J = 1, 4 ) / 2573, 1368, 3508, \$ 3017 / DATA ( MM( 127, J ), J = 1, 4 ) / 1148, 1848, 3525, \$ 2141 / DATA ( MM( 128, J ), J = 1, 4 ) / 545, 2366, 3801, \$ 1537 / * .. * .. Executable Statements .. * I1 = ISEED( 1 ) I2 = ISEED( 2 ) I3 = ISEED( 3 ) I4 = ISEED( 4 ) * DO 10 I = 1, MIN( N, LV ) * * Multiply the seed by i-th power of the multiplier modulo 2**48 * IT4 = I4*MM( I, 4 ) IT3 = IT4 / IPW2 IT4 = IT4 - IPW2*IT3 IT3 = IT3 + I3*MM( I, 4 ) + I4*MM( I, 3 ) IT2 = IT3 / IPW2 IT3 = IT3 - IPW2*IT2 IT2 = IT2 + I2*MM( I, 4 ) + I3*MM( I, 3 ) + I4*MM( I, 2 ) IT1 = IT2 / IPW2 IT2 = IT2 - IPW2*IT1 IT1 = IT1 + I1*MM( I, 4 ) + I2*MM( I, 3 ) + I3*MM( I, 2 ) + \$ I4*MM( I, 1 ) IT1 = MOD( IT1, IPW2 ) * * Convert 48-bit integer to a real number in the interval (0,1) * X( I ) = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* \$ DBLE( IT4 ) ) ) ) 10 CONTINUE * * Return final value of seed * ISEED( 1 ) = IT1 ISEED( 2 ) = IT2 ISEED( 3 ) = IT3 ISEED( 4 ) = IT4 RETURN * * End of DLARUV * END