![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DLARAN 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 * Definition: 00009 * =========== 00010 * 00011 * DOUBLE PRECISION FUNCTION DLARAN( ISEED ) 00012 * 00013 * .. Array Arguments .. 00014 * INTEGER ISEED( 4 ) 00015 * .. 00016 * 00017 * 00018 *> \par Purpose: 00019 * ============= 00020 *> 00021 *> \verbatim 00022 *> 00023 *> DLARAN returns a random real number from a uniform (0,1) 00024 *> distribution. 00025 *> \endverbatim 00026 * 00027 * Arguments: 00028 * ========== 00029 * 00030 *> \param[in,out] ISEED 00031 *> \verbatim 00032 *> ISEED is INTEGER array, dimension (4) 00033 *> On entry, the seed of the random number generator; the array 00034 *> elements must be between 0 and 4095, and ISEED(4) must be 00035 *> odd. 00036 *> On exit, the seed is updated. 00037 *> \endverbatim 00038 * 00039 * Authors: 00040 * ======== 00041 * 00042 *> \author Univ. of Tennessee 00043 *> \author Univ. of California Berkeley 00044 *> \author Univ. of Colorado Denver 00045 *> \author NAG Ltd. 00046 * 00047 *> \date November 2011 00048 * 00049 *> \ingroup list_temp 00050 * 00051 *> \par Further Details: 00052 * ===================== 00053 *> 00054 *> \verbatim 00055 *> 00056 *> This routine uses a multiplicative congruential method with modulus 00057 *> 2**48 and multiplier 33952834046453 (see G.S.Fishman, 00058 *> 'Multiplicative congruential random number generators with modulus 00059 *> 2**b: an exhaustive analysis for b = 32 and a partial analysis for 00060 *> b = 48', Math. Comp. 189, pp 331-344, 1990). 00061 *> 00062 *> 48-bit integers are stored in 4 integer array elements with 12 bits 00063 *> per element. Hence the routine is portable across machines with 00064 *> integers of 32 bits or more. 00065 *> \endverbatim 00066 *> 00067 * ===================================================================== 00068 DOUBLE PRECISION FUNCTION DLARAN( ISEED ) 00069 * 00070 * -- LAPACK auxiliary routine (version 3.4.0) -- 00071 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00072 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00073 * November 2011 00074 * 00075 * .. Array Arguments .. 00076 INTEGER ISEED( 4 ) 00077 * .. 00078 * 00079 * ===================================================================== 00080 * 00081 * .. Parameters .. 00082 INTEGER M1, M2, M3, M4 00083 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 ) 00084 DOUBLE PRECISION ONE 00085 PARAMETER ( ONE = 1.0D+0 ) 00086 INTEGER IPW2 00087 DOUBLE PRECISION R 00088 PARAMETER ( IPW2 = 4096, R = ONE / IPW2 ) 00089 * .. 00090 * .. Local Scalars .. 00091 INTEGER IT1, IT2, IT3, IT4 00092 DOUBLE PRECISION RNDOUT 00093 * .. 00094 * .. Intrinsic Functions .. 00095 INTRINSIC DBLE, MOD 00096 * .. 00097 * .. Executable Statements .. 00098 10 CONTINUE 00099 * 00100 * multiply the seed by the multiplier modulo 2**48 00101 * 00102 IT4 = ISEED( 4 )*M4 00103 IT3 = IT4 / IPW2 00104 IT4 = IT4 - IPW2*IT3 00105 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3 00106 IT2 = IT3 / IPW2 00107 IT3 = IT3 - IPW2*IT2 00108 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2 00109 IT1 = IT2 / IPW2 00110 IT2 = IT2 - IPW2*IT1 00111 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 + 00112 $ ISEED( 4 )*M1 00113 IT1 = MOD( IT1, IPW2 ) 00114 * 00115 * return updated seed 00116 * 00117 ISEED( 1 ) = IT1 00118 ISEED( 2 ) = IT2 00119 ISEED( 3 ) = IT3 00120 ISEED( 4 ) = IT4 00121 * 00122 * convert 48-bit integer to a real number in the interval (0,1) 00123 * 00124 RNDOUT = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* 00125 $ ( DBLE( IT4 ) ) ) ) ) 00126 * 00127 IF (RNDOUT.EQ.1.0D+0) THEN 00128 * If a real number has n bits of precision, and the first 00129 * n bits of the 48-bit integer above happen to be all 1 (which 00130 * will occur about once every 2**n calls), then DLARAN will 00131 * be rounded to exactly 1.0. 00132 * Since DLARAN is not supposed to return exactly 0.0 or 1.0 00133 * (and some callers of DLARAN, such as CLARND, depend on that), 00134 * the statistically correct thing to do in this situation is 00135 * simply to iterate again. 00136 * N.B. the case DLARAN = 0.0 should not be possible. 00137 * 00138 GOTO 10 00139 END IF 00140 * 00141 DLARAN = RNDOUT 00142 RETURN 00143 * 00144 * End of DLARAN 00145 * 00146 END