Télécharger dlaruv.eso

Retour à la liste

Numérotation des lignes :

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

© Cast3M 2003 - Tous droits réservés.
Mentions légales