SRC/dlaruv.f

   1 *> \brief \b DLARUV 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 DLARUV + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dlaruv.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dlaruv.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dlaruv.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DLARUV( ISEED, N, X )
  22 *
  23 *       .. Scalar Arguments ..
  24 *       INTEGER            N
  25 *       ..
  26 *       .. Array Arguments ..
  27 *       INTEGER            ISEED( 4 )
  28 *       DOUBLE PRECISION   X( N )
  29 *       ..
  30 *
  31 *
  32 *> \par Purpose:
  33 *  =============
  34 *>
  35 *> \verbatim
  36 *>
  37 *> DLARUV 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 DLARNV and ZLARNV.
  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 DOUBLE PRECISION 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 OTHERauxiliary
  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 DLARUV( 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       DOUBLE PRECISION   X( N )
 109 *     ..
 110 *
 111 *  =====================================================================
 112 *
 113 *     .. Parameters ..
 114       DOUBLE PRECISION   ONE
 115       PARAMETER          ( ONE = 1.0D0 )
 116       INTEGER            LV, IPW2
 117       DOUBLE PRECISION   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          DBLE, MIN, MOD
 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*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
 416      $            DBLE( IT4 ) ) ) )
 417 *
 418          IF (X( I ).EQ.1.0D0) 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.
 423 *           Since X( I ) is not supposed to return exactly 0.0 or 1.0,
 424 *           the statistically correct thing to do in this situation is
 425 *           simply to iterate again.
 426 *           N.B. the case X( I ) = 0.0 should not be possible.
 427             I1 = I1 + 2
 428             I2 = I2 + 2
 429             I3 = I3 + 2
 430             I4 = I4 + 2
 431             GOTO 20
 432          END IF
 433 *
 434    10 CONTINUE
 435 *
 436 *     Return final value of seed
 437 *
 438       ISEED( 1 ) = IT1
 439       ISEED( 2 ) = IT2
 440       ISEED( 3 ) = IT3
 441       ISEED( 4 ) = IT4
 442       RETURN
 443 *
 444 *     End of DLARUV
 445 *
 446       END