*> \brief \b DLARAN * * =========== DOCUMENTATION =========== * * Online html documentation available at * http://www.netlib.org/lapack/explore-html/ * * Definition: * =========== * * DOUBLE PRECISION FUNCTION DLARAN( ISEED ) * * .. Array Arguments .. * INTEGER ISEED( 4 ) * .. * * *> \par Purpose: * ============= *> *> \verbatim *> *> DLARAN returns a random real number from a uniform (0,1) *> distribution. *> \endverbatim * * Arguments: * ========== * *> \param[in,out] ISEED *> \verbatim *> ISEED is INTEGER array, dimension (4) *> On entry, the seed of the random number generator; the array *> elements must be between 0 and 4095, and ISEED(4) must be *> odd. *> On exit, the seed is updated. *> \endverbatim * * Authors: * ======== * *> \author Univ. of Tennessee *> \author Univ. of California Berkeley *> \author Univ. of Colorado Denver *> \author NAG Ltd. * *> \date November 2011 * *> \ingroup list_temp * *> \par Further Details: * ===================== *> *> \verbatim *> *> This routine uses a multiplicative congruential method with modulus *> 2**48 and multiplier 33952834046453 (see G.S.Fishman, *> 'Multiplicative congruential random number generators with modulus *> 2**b: an exhaustive analysis for b = 32 and a partial analysis for *> b = 48', Math. Comp. 189, pp 331-344, 1990). *> *> 48-bit integers are stored in 4 integer array elements with 12 bits *> per element. Hence the routine is portable across machines with *> integers of 32 bits or more. *> \endverbatim *> * ===================================================================== DOUBLE PRECISION FUNCTION DLARAN( ISEED ) * * -- LAPACK auxiliary routine (version 3.4.0) -- * -- LAPACK is a software package provided by Univ. of Tennessee, -- * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- * November 2011 * * .. Array Arguments .. INTEGER ISEED( 4 ) * .. * * ===================================================================== * * .. Parameters .. INTEGER M1, M2, M3, M4 PARAMETER ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 ) DOUBLE PRECISION ONE PARAMETER ( ONE = 1.0D+0 ) INTEGER IPW2 DOUBLE PRECISION R PARAMETER ( IPW2 = 4096, R = ONE / IPW2 ) * .. * .. Local Scalars .. INTEGER IT1, IT2, IT3, IT4 DOUBLE PRECISION RNDOUT * .. * .. Intrinsic Functions .. INTRINSIC DBLE, MOD * .. * .. Executable Statements .. 10 CONTINUE * * multiply the seed by the multiplier modulo 2**48 * IT4 = ISEED( 4 )*M4 IT3 = IT4 / IPW2 IT4 = IT4 - IPW2*IT3 IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3 IT2 = IT3 / IPW2 IT3 = IT3 - IPW2*IT2 IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2 IT1 = IT2 / IPW2 IT2 = IT2 - IPW2*IT1 IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 + $ ISEED( 4 )*M1 IT1 = MOD( IT1, IPW2 ) * * return updated seed * ISEED( 1 ) = IT1 ISEED( 2 ) = IT2 ISEED( 3 ) = IT3 ISEED( 4 ) = IT4 * * convert 48-bit integer to a real number in the interval (0,1) * RNDOUT = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R* $ ( DBLE( IT4 ) ) ) ) ) * IF (RNDOUT.EQ.1.0D+0) THEN * If a real number has n bits of precision, and the first * n bits of the 48-bit integer above happen to be all 1 (which * will occur about once every 2**n calls), then DLARAN will * be rounded to exactly 1.0. * Since DLARAN is not supposed to return exactly 0.0 or 1.0 * (and some callers of DLARAN, such as CLARND, depend on that), * the statistically correct thing to do in this situation is * simply to iterate again. * N.B. the case DLARAN = 0.0 should not be possible. * GOTO 10 END IF * DLARAN = RNDOUT RETURN * * End of DLARAN * END