kusano 7d535a
#include "f2c.h"
kusano 7d535a
kusano 7d535a
/* Subroutine */ int zlarnv_(integer *idist, integer *iseed, integer *n, 
kusano 7d535a
	doublecomplex *x)
kusano 7d535a
{
kusano 7d535a
/*  -- LAPACK auxiliary routine (version 2.0) --   
kusano 7d535a
       Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,   
kusano 7d535a
       Courant Institute, Argonne National Lab, and Rice University   
kusano 7d535a
       September 30, 1994   
kusano 7d535a
kusano 7d535a
kusano 7d535a
    Purpose   
kusano 7d535a
    =======   
kusano 7d535a
kusano 7d535a
    ZLARNV returns a vector of n random complex numbers from a uniform or 
kusano 7d535a
  
kusano 7d535a
    normal distribution.   
kusano 7d535a
kusano 7d535a
    Arguments   
kusano 7d535a
    =========   
kusano 7d535a
kusano 7d535a
    IDIST   (input) INTEGER   
kusano 7d535a
            Specifies the distribution of the random numbers:   
kusano 7d535a
            = 1:  real and imaginary parts each uniform (0,1)   
kusano 7d535a
            = 2:  real and imaginary parts each uniform (-1,1)   
kusano 7d535a
            = 3:  real and imaginary parts each normal (0,1)   
kusano 7d535a
            = 4:  uniformly distributed on the disc abs(z) < 1   
kusano 7d535a
            = 5:  uniformly distributed on the circle abs(z) = 1   
kusano 7d535a
kusano 7d535a
    ISEED   (input/output) INTEGER array, dimension (4)   
kusano 7d535a
            On entry, the seed of the random number generator; the array 
kusano 7d535a
  
kusano 7d535a
            elements must be between 0 and 4095, and ISEED(4) must be   
kusano 7d535a
            odd.   
kusano 7d535a
            On exit, the seed is updated.   
kusano 7d535a
kusano 7d535a
    N       (input) INTEGER   
kusano 7d535a
            The number of random numbers to be generated.   
kusano 7d535a
kusano 7d535a
    X       (output) COMPLEX*16 array, dimension (N)   
kusano 7d535a
            The generated random numbers.   
kusano 7d535a
kusano 7d535a
    Further Details   
kusano 7d535a
    ===============   
kusano 7d535a
kusano 7d535a
    This routine calls the auxiliary routine DLARUV to generate random   
kusano 7d535a
    real numbers from a uniform (0,1) distribution, in batches of up to   
kusano 7d535a
    128 using vectorisable code. The Box-Muller method is used to   
kusano 7d535a
    transform numbers from a uniform to a normal distribution.   
kusano 7d535a
kusano 7d535a
    ===================================================================== 
kusano 7d535a
  
kusano 7d535a
kusano 7d535a
kusano 7d535a
    
kusano 7d535a
   Parameter adjustments   
kusano 7d535a
       Function Body */
kusano 7d535a
    /* System generated locals */
kusano 7d535a
    integer i__1, i__2, i__3, i__4, i__5;
kusano 7d535a
    doublereal d__1, d__2;
kusano 7d535a
    doublecomplex z__1, z__2, z__3;
kusano 7d535a
    /* Builtin functions */
kusano 7d535a
    double log(doublereal), sqrt(doublereal);
kusano 7d535a
    void z_exp(doublecomplex *, doublecomplex *);
kusano 7d535a
    /* Local variables */
kusano 7d535a
    static integer i;
kusano 7d535a
    static doublereal u[128];
kusano 7d535a
    static integer il, iv;
kusano 7d535a
    extern /* Subroutine */ int dlaruv_(integer *, integer *, doublereal *);
kusano 7d535a
kusano 7d535a
kusano 7d535a
#define U(I) u[(I)]
kusano 7d535a
#define X(I) x[(I)-1]
kusano 7d535a
#define ISEED(I) iseed[(I)-1]
kusano 7d535a
kusano 7d535a
kusano 7d535a
    i__1 = *n;
kusano 7d535a
    for (iv = 1; iv <= *n; iv += 64) {
kusano 7d535a
/* Computing MIN */
kusano 7d535a
	i__2 = 64, i__3 = *n - iv + 1;
kusano 7d535a
	il = min(i__2,i__3);
kusano 7d535a
kusano 7d535a
/*        Call DLARUV to generate 2*IL real numbers from a uniform (0,
kusano 7d535a
1)   
kusano 7d535a
          distribution (2*IL <= LV) */
kusano 7d535a
kusano 7d535a
	i__2 = il << 1;
kusano 7d535a
	dlaruv_(&ISEED(1), &i__2, u);
kusano 7d535a
kusano 7d535a
	if (*idist == 1) {
kusano 7d535a
kusano 7d535a
/*           Copy generated numbers */
kusano 7d535a
kusano 7d535a
	    i__2 = il;
kusano 7d535a
	    for (i = 1; i <= il; ++i) {
kusano 7d535a
		i__3 = iv + i - 1;
kusano 7d535a
		i__4 = (i << 1) - 2;
kusano 7d535a
		i__5 = (i << 1) - 1;
kusano 7d535a
		z__1.r = U((i<<1)-2), z__1.i = U((i<<1)-1);
kusano 7d535a
		X(iv+i-1).r = z__1.r, X(iv+i-1).i = z__1.i;
kusano 7d535a
/* L10: */
kusano 7d535a
	    }
kusano 7d535a
	} else if (*idist == 2) {
kusano 7d535a
kusano 7d535a
/*           Convert generated numbers to uniform (-1,1) distribut
kusano 7d535a
ion */
kusano 7d535a
kusano 7d535a
	    i__2 = il;
kusano 7d535a
	    for (i = 1; i <= il; ++i) {
kusano 7d535a
		i__3 = iv + i - 1;
kusano 7d535a
		d__1 = U((i << 1) - 2) * 2. - 1.;
kusano 7d535a
		d__2 = U((i << 1) - 1) * 2. - 1.;
kusano 7d535a
		z__1.r = d__1, z__1.i = d__2;
kusano 7d535a
		X(iv+i-1).r = z__1.r, X(iv+i-1).i = z__1.i;
kusano 7d535a
/* L20: */
kusano 7d535a
	    }
kusano 7d535a
	} else if (*idist == 3) {
kusano 7d535a
kusano 7d535a
/*           Convert generated numbers to normal (0,1) distributio
kusano 7d535a
n */
kusano 7d535a
kusano 7d535a
	    i__2 = il;
kusano 7d535a
	    for (i = 1; i <= il; ++i) {
kusano 7d535a
		i__3 = iv + i - 1;
kusano 7d535a
		d__1 = sqrt(log(U((i << 1) - 2)) * -2.);
kusano 7d535a
		d__2 = U((i << 1) - 1) * 6.2831853071795864769252867663;
kusano 7d535a
		z__3.r = 0., z__3.i = d__2;
kusano 7d535a
		z_exp(&z__2, &z__3);
kusano 7d535a
		z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
kusano 7d535a
		X(iv+i-1).r = z__1.r, X(iv+i-1).i = z__1.i;
kusano 7d535a
/* L30: */
kusano 7d535a
	    }
kusano 7d535a
	} else if (*idist == 4) {
kusano 7d535a
kusano 7d535a
/*           Convert generated numbers to complex numbers uniforml
kusano 7d535a
y   
kusano 7d535a
             distributed on the unit disk */
kusano 7d535a
kusano 7d535a
	    i__2 = il;
kusano 7d535a
	    for (i = 1; i <= il; ++i) {
kusano 7d535a
		i__3 = iv + i - 1;
kusano 7d535a
		d__1 = sqrt(U((i << 1) - 2));
kusano 7d535a
		d__2 = U((i << 1) - 1) * 6.2831853071795864769252867663;
kusano 7d535a
		z__3.r = 0., z__3.i = d__2;
kusano 7d535a
		z_exp(&z__2, &z__3);
kusano 7d535a
		z__1.r = d__1 * z__2.r, z__1.i = d__1 * z__2.i;
kusano 7d535a
		X(iv+i-1).r = z__1.r, X(iv+i-1).i = z__1.i;
kusano 7d535a
/* L40: */
kusano 7d535a
	    }
kusano 7d535a
	} else if (*idist == 5) {
kusano 7d535a
kusano 7d535a
/*           Convert generated numbers to complex numbers uniforml
kusano 7d535a
y   
kusano 7d535a
             distributed on the unit circle */
kusano 7d535a
kusano 7d535a
	    i__2 = il;
kusano 7d535a
	    for (i = 1; i <= il; ++i) {
kusano 7d535a
		i__3 = iv + i - 1;
kusano 7d535a
		d__1 = U((i << 1) - 1) * 6.2831853071795864769252867663;
kusano 7d535a
		z__2.r = 0., z__2.i = d__1;
kusano 7d535a
		z_exp(&z__1, &z__2);
kusano 7d535a
		X(iv+i-1).r = z__1.r, X(iv+i-1).i = z__1.i;
kusano 7d535a
/* L50: */
kusano 7d535a
	    }
kusano 7d535a
	}
kusano 7d535a
/* L60: */
kusano 7d535a
    }
kusano 7d535a
    return 0;
kusano 7d535a
kusano 7d535a
/*     End of ZLARNV */
kusano 7d535a
kusano 7d535a
} /* zlarnv_ */
kusano 7d535a