LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlaran.f
Go to the documentation of this file.
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
 All Files Functions