git.maemo.org Git - opencv/blob - 3rdparty/lapack/slarnv.c

   1 #include "clapack.h"
   2
   3 /* Subroutine */ int slarnv_(integer *idist, integer *iseed, integer *n, real
   4         *x)
   5 {
   6     /* System generated locals */
   7     integer i__1, i__2, i__3;
   8
   9     /* Builtin functions */
  10     double log(doublereal), sqrt(doublereal), cos(doublereal);
  11
  12     /* Local variables */
  13     integer i__;
  14     real u[128];
  15     integer il, iv, il2;
  16     extern /* Subroutine */ int slaruv_(integer *, integer *, real *);
  17
  18
  19 /*  -- LAPACK auxiliary routine (version 3.1) -- */
  20 /*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. */
  21 /*     November 2006 */
  22
  23 /*     .. Scalar Arguments .. */
  24 /*     .. */
  25 /*     .. Array Arguments .. */
  26 /*     .. */
  27
  28 /*  Purpose */
  29 /*  ======= */
  30
  31 /*  SLARNV returns a vector of n random real numbers from a uniform or */
  32 /*  normal distribution. */
  33
  34 /*  Arguments */
  35 /*  ========= */
  36
  37 /*  IDIST   (input) INTEGER */
  38 /*          Specifies the distribution of the random numbers: */
  39 /*          = 1:  uniform (0,1) */
  40 /*          = 2:  uniform (-1,1) */
  41 /*          = 3:  normal (0,1) */
  42
  43 /*  ISEED   (input/output) INTEGER array, dimension (4) */
  44 /*          On entry, the seed of the random number generator; the array */
  45 /*          elements must be between 0 and 4095, and ISEED(4) must be */
  46 /*          odd. */
  47 /*          On exit, the seed is updated. */
  48
  49 /*  N       (input) INTEGER */
  50 /*          The number of random numbers to be generated. */
  51
  52 /*  X       (output) REAL array, dimension (N) */
  53 /*          The generated random numbers. */
  54
  55 /*  Further Details */
  56 /*  =============== */
  57
  58 /*  This routine calls the auxiliary routine SLARUV to generate random */
  59 /*  real numbers from a uniform (0,1) distribution, in batches of up to */
  60 /*  128 using vectorisable code. The Box-Muller method is used to */
  61 /*  transform numbers from a uniform to a normal distribution. */
  62
  63 /*  ===================================================================== */
  64
  65 /*     .. Parameters .. */
  66 /*     .. */
  67 /*     .. Local Scalars .. */
  68 /*     .. */
  69 /*     .. Local Arrays .. */
  70 /*     .. */
  71 /*     .. Intrinsic Functions .. */
  72 /*     .. */
  73 /*     .. External Subroutines .. */
  74 /*     .. */
  75 /*     .. Executable Statements .. */
  76
  77     /* Parameter adjustments */
  78     --x;
  79     --iseed;
  80
  81     /* Function Body */
  82     i__1 = *n;
  83     for (iv = 1; iv <= i__1; iv += 64) {
  84 /* Computing MIN */
  85         i__2 = 64, i__3 = *n - iv + 1;
  86         il = min(i__2,i__3);
  87         if (*idist == 3) {
  88             il2 = il << 1;
  89         } else {
  90             il2 = il;
  91         }
  92
  93 /*        Call SLARUV to generate IL2 numbers from a uniform (0,1) */
  94 /*        distribution (IL2 <= LV) */
  95
  96         slaruv_(&iseed[1], &il2, u);
  97
  98         if (*idist == 1) {
  99
 100 /*           Copy generated numbers */
 101
 102             i__2 = il;
 103             for (i__ = 1; i__ <= i__2; ++i__) {
 104                 x[iv + i__ - 1] = u[i__ - 1];
 105 /* L10: */
 106             }
 107         } else if (*idist == 2) {
 108
 109 /*           Convert generated numbers to uniform (-1,1) distribution */
 110
 111             i__2 = il;
 112             for (i__ = 1; i__ <= i__2; ++i__) {
 113                 x[iv + i__ - 1] = u[i__ - 1] * 2.f - 1.f;
 114 /* L20: */
 115             }
 116         } else if (*idist == 3) {
 117
 118 /*           Convert generated numbers to normal (0,1) distribution */
 119
 120             i__2 = il;
 121             for (i__ = 1; i__ <= i__2; ++i__) {
 122                 x[iv + i__ - 1] = sqrt(log(u[(i__ << 1) - 2]) * -2.f) * cos(u[
 123                         (i__ << 1) - 1] * 6.2831853071795864769252867663f);
 124 /* L30: */
 125             }
 126         }
 127 /* L40: */
 128     }
 129     return 0;
 130
 131 /*     End of SLARNV */
 132
 133 } /* slarnv_ */