Return to dlaruv.f CVS log

Up to [local] / rpl / lapack / lapack

Cohérence.

    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 auxOTHERauxiliary
   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