LAPACK 3.3.1
Linear Algebra PACKage

slaruv.f

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