TESTING/MATGEN/dlaran.f

   1 *> \brief \b DLARAN
   2 *
   3 *  =========== DOCUMENTATION ===========
   4 *
   5 * Online html documentation available at
   6 *            http://www.netlib.org/lapack/explore-html/
   7 *
   8 *  Definition:
   9 *  ===========
  10 *
  11 *       DOUBLE PRECISION FUNCTION DLARAN( ISEED )
  12 *
  13 *       .. Array Arguments ..
  14 *       INTEGER            ISEED( 4 )
  15 *       ..
  16 *
  17 *
  18 *> \par Purpose:
  19 *  =============
  20 *>
  21 *> \verbatim
  22 *>
  23 *> DLARAN returns a random real number from a uniform (0,1)
  24 *> distribution.
  25 *> \endverbatim
  26 *
  27 *  Arguments:
  28 *  ==========
  29 *
  30 *> \param[in,out] ISEED
  31 *> \verbatim
  32 *>          ISEED is INTEGER array, dimension (4)
  33 *>          On entry, the seed of the random number generator; the array
  34 *>          elements must be between 0 and 4095, and ISEED(4) must be
  35 *>          odd.
  36 *>          On exit, the seed is updated.
  37 *> \endverbatim
  38 *
  39 *  Authors:
  40 *  ========
  41 *
  42 *> \author Univ. of Tennessee
  43 *> \author Univ. of California Berkeley
  44 *> \author Univ. of Colorado Denver
  45 *> \author NAG Ltd.
  46 *
  47 *> \date November 2011
  48 *
  49 *> \ingroup list_temp
  50 *
  51 *> \par Further Details:
  52 *  =====================
  53 *>
  54 *> \verbatim
  55 *>
  56 *>  This routine uses a multiplicative congruential method with modulus
  57 *>  2**48 and multiplier 33952834046453 (see G.S.Fishman,
  58 *>  'Multiplicative congruential random number generators with modulus
  59 *>  2**b: an exhaustive analysis for b = 32 and a partial analysis for
  60 *>  b = 48', Math. Comp. 189, pp 331-344, 1990).
  61 *>
  62 *>  48-bit integers are stored in 4 integer array elements with 12 bits
  63 *>  per element. Hence the routine is portable across machines with
  64 *>  integers of 32 bits or more.
  65 *> \endverbatim
  66 *>
  67 *  =====================================================================
  68       DOUBLE PRECISION FUNCTION DLARAN( ISEED )
  69 *
  70 *  -- LAPACK auxiliary routine (version 3.4.0) --
  71 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
  72 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
  73 *     November 2011
  74 *
  75 *     .. Array Arguments ..
  76       INTEGER            ISEED( 4 )
  77 *     ..
  78 *
  79 *  =====================================================================
  80 *
  81 *     .. Parameters ..
  82       INTEGER            M1, M2, M3, M4
  83       PARAMETER          ( M1 = 494, M2 = 322, M3 = 2508, M4 = 2549 )
  84       DOUBLE PRECISION   ONE
  85       PARAMETER          ( ONE = 1.0D+0 )
  86       INTEGER            IPW2
  87       DOUBLE PRECISION   R
  88       PARAMETER          ( IPW2 = 4096, R = ONE / IPW2 )
  89 *     ..
  90 *     .. Local Scalars ..
  91       INTEGER            IT1, IT2, IT3, IT4
  92       DOUBLE PRECISION   RNDOUT
  93 *     ..
  94 *     .. Intrinsic Functions ..
  95       INTRINSIC          DBLE, MOD
  96 *     ..
  97 *     .. Executable Statements ..
  98   10  CONTINUE
  99 *
 100 *     multiply the seed by the multiplier modulo 2**48
 101 *
 102       IT4 = ISEED( 4 )*M4
 103       IT3 = IT4 / IPW2
 104       IT4 = IT4 - IPW2*IT3
 105       IT3 = IT3 + ISEED( 3 )*M4 + ISEED( 4 )*M3
 106       IT2 = IT3 / IPW2
 107       IT3 = IT3 - IPW2*IT2
 108       IT2 = IT2 + ISEED( 2 )*M4 + ISEED( 3 )*M3 + ISEED( 4 )*M2
 109       IT1 = IT2 / IPW2
 110       IT2 = IT2 - IPW2*IT1
 111       IT1 = IT1 + ISEED( 1 )*M4 + ISEED( 2 )*M3 + ISEED( 3 )*M2 +
 112      $      ISEED( 4 )*M1
 113       IT1 = MOD( IT1, IPW2 )
 114 *
 115 *     return updated seed
 116 *
 117       ISEED( 1 ) = IT1
 118       ISEED( 2 ) = IT2
 119       ISEED( 3 ) = IT3
 120       ISEED( 4 ) = IT4
 121 *
 122 *     convert 48-bit integer to a real number in the interval (0,1)
 123 *
 124       RNDOUT = R*( DBLE( IT1 )+R*( DBLE( IT2 )+R*( DBLE( IT3 )+R*
 125      $         ( DBLE( IT4 ) ) ) ) )
 126 *
 127       IF (RNDOUT.EQ.1.0D+0) THEN
 128 *        If a real number has n bits of precision, and the first
 129 *        n bits of the 48-bit integer above happen to be all 1 (which
 130 *        will occur about once every 2**n calls), then DLARAN will
 131 *        be rounded to exactly 1.0.
 132 *        Since DLARAN is not supposed to return exactly 0.0 or 1.0
 133 *        (and some callers of DLARAN, such as CLARND, depend on that),
 134 *        the statistically correct thing to do in this situation is
 135 *        simply to iterate again.
 136 *        N.B. the case DLARAN = 0.0 should not be possible.
 137 *
 138          GOTO 10
 139       END IF
 140 *
 141       DLARAN = RNDOUT
 142       RETURN
 143 *
 144 *     End of DLARAN
 145 *
 146       END