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