Header for random number configuration. More...

#include <math.h>
#include <par_mpi.h>
#include <sizes.h>

Include dependency graph for random.h:

This graph shows which files directly or indirectly include this file:

Macros
#define	RAN2

#define	IM1 2147483563

#define	IM2 2147483399

#define	AM (1.0/IM1)

#define	IMM1 (IM1-1)

#define	IA1 40014

#define	IA2 40692

#define	IQ1 53668

#define	IQ2 52774

#define	IR1 12211

#define	IR2 3791

#define	NTAB 32

#define	NDIV (1+IMM1/NTAB)

#define	EPS 1.2e-7

#define	RNMX (1.0-EPS)

Functions
int	ranset (long *seed)
	Dummy seed the ran2 generator.

int	Par_ranset (long *seed, int iread)
	Uses the rank to get a new seed. Copying from the FORTRAN description here c create new seeds in range seed to 9seed c having a range of 0seed gave an unfortunate pattern c in the underlying value of ds(1) (it was always 10 times bigger c on the last processor). This does not appear to happen with 9.

double	ran2 (long *idum)
	Generates uniformly distributed random double between zero and one as described in numerical recipes. It's also thread-safe for different seeds.

int	Gauss_z (Complex *ps, unsigned int n, const Complex mu, const double sigma)
	Generates a vector of normally distributed random double precision complex numbers using the Box-Muller Method.

int	Gauss_d (double *ps, unsigned int n, const double mu, const double sigma)
	Generates a vector of normally distributed random double precision numbers using the Box-Muller Method.

int	Gauss_c (Complex_f *ps, unsigned int n, const Complex_f mu, const float sigma)
	Generates a vector of normally distributed random single precision complex numbers using the Box-Muller Method.

int	Gauss_f (float *ps, unsigned int n, const float mu, const float sigma)
	Generates a vector of normally distributed random single precision numbers using the Box-Muller Method.

int	Par_ranread (char filename, double ranval)
	Reads ps from a file Since this function is very similar to Par_sread, I'm not really going to comment it check there if you are confused about things.

double	Par_granf ()
	Generates a random double which is then sent to the other ranks.

int	ran_test ()
	Test Functions.

Variables
long	seed
	RAN2 seed.

Detailed Description

Header for random number configuration.

Definition in file random.h.

Macro Definition Documentation

◆ AM

#define AM (1.0/IM1)

Definition at line 234 of file random.h.

◆ EPS

#define EPS 1.2e-7

Definition at line 244 of file random.h.

◆ IA1

#define IA1 40014

Definition at line 236 of file random.h.

◆ IA2

#define IA2 40692

Definition at line 237 of file random.h.

◆ IM1

#define IM1 2147483563

Definition at line 232 of file random.h.

◆ IM2

#define IM2 2147483399

Definition at line 233 of file random.h.

◆ IMM1

#define IMM1 (IM1-1)

Definition at line 235 of file random.h.

◆ IQ1

#define IQ1 53668

Definition at line 238 of file random.h.

◆ IQ2

#define IQ2 52774

Definition at line 239 of file random.h.

◆ IR1

#define IR1 12211

Definition at line 240 of file random.h.

◆ IR2

#define IR2 3791

Definition at line 241 of file random.h.

◆ NDIV

#define NDIV (1+IMM1/NTAB)

Definition at line 243 of file random.h.

◆ NTAB

#define NTAB 32

Definition at line 242 of file random.h.

◆ RAN2

#define RAN2

Definition at line 231 of file random.h.

◆ RNMX

#define RNMX (1.0-EPS)

Definition at line 245 of file random.h.

Function Documentation

◆ Gauss_c()

int Gauss_c	(	Complex_f *	ps,
		unsigned int	n,
		const Complex_f	mu,
		const float	sigma )

Generates a vector of normally distributed random single precision complex numbers using the Box-Muller Method.

Parameters

ps	The output array
n	The array length
mu	mean
sigma	variance

Returns: Zero on success integer error code otherwise

Definition at line 260 of file random.c.

                                                                                 {
   /*
    * @brief   Generates a vector of normally distributed random single precision complex numbers using the Box-Muller Method
    * 
    * @param   ps:      The output array
    * @param   n:       The array length
    * @param   mu:      mean
    * @param   sigma:   variance
    *
    * @see ran2(), par_dcopy(), gsl_rng_uniform(), vdRngUniform()
    * 
    * @return Zero on success integer error code otherwise
    */
   const char *funcname = "Gauss_z";
   if(n<=0){
      fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
            ARRAYLEN, funcname, n);
#if(nproc>1)
      MPI_Abort(comm,ARRAYLEN);
#else
      exit(ARRAYLEN);
#endif
   }
#pragma unroll
   for(int i=0;i<n;i++){
      /* Marsaglia Method for fun
         do{
         u=sfmt_genrand_real1(sfmt);
         v=sfmt_genrand_real1(sfmt);
         r=u*u+v*v;
         }while(0<r & r<1);
         r=sqrt(r);
         r=sqrt(-2.0*log(r)/r)*sigma;
         ps[i] = mu+u*r + I*(mu+v*r);
       */
#ifdef __RANLUX__
      float r =sigma*sqrt(-2*log(gsl_rng_uniform(ranlux_instd)));
      float theta=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
#else
      float r =sigma*sqrt(-2*log(ran2(&seed)));
      float theta=2.0*M_PI*ran2(&seed);
#endif
      ps[i]=r*(cos(theta)+mu+sin(theta)*I)+mu;
   }     
   return 0;
}

References M_PI, ran2(), and seed.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Gauss_d()

int Gauss_d	(	double *	ps,
		unsigned int	n,
		const double	mu,
		const double	sigma )

Generates a vector of normally distributed random double precision numbers using the Box-Muller Method.

Parameters

ps	The output array
n	The array length
mu	mean
sigma	variance

Returns: Zero on success integer error code otherwise

Definition at line 306 of file random.c.

                                                                            {
   /*
    * @brief   Generates a vector of normally distributed random double precision numbers using the Box-Muller Method
    * 
    * @param   ps:      The output array
    * @param   n:       The array length
    * @param   mu:      mean
    * @param   sigma:   variance
    *
    * @see ran2(), par_dcopy(), gsl_rng_uniform(), vdRngUniform()
    * 
    * @return Zero on success integer error code otherwise
    */
   const char *funcname = "Gauss_z";
   //The FORTRAN Code had two different Gauss Routines. gaussp having unit
   //mean and variance and gauss0 where the variance would appear to be 1/sqrt(2)
   //(Since we multiply by sqrt(-ln(r)) instead of sqrt(-2ln(r)) )
   if(n<=0){
      fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
            ARRAYLEN, funcname, n);
#if(nproc>1)
      MPI_Abort(comm,ARRAYLEN);
#else
      exit(ARRAYLEN);
#endif
   }
   int i;
   double r, u, v;
   //If n is odd we calculate the last index seperately and the rest in pairs
   if(n%2==1){
      n--;
#ifdef __RANLUX__
      r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
      ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
#else
      r=2.0*M_PI*ran2(&seed);
      ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
#endif
   }
   for(i=0;i<n;i+=2){
      /* Marsaglia Method for fun
         do{
         u=sfmt_genrand_real1(sfmt);
         v=sfmt_genrand_real1(sfmt);
         r=u*u+v*v;
         }while(0<r & r<1);
         r=sqrt(r);
         r=sqrt(-2.0*log(r)/r)*sigma;
         ps[i] = mu+u*r; 
         ps[i+1]=mu+v*r;
       */
#ifdef __RANLUX__
      u=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*sigma;
      r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
#else
      u=sqrt(-2*log(ran2(&seed)))*sigma;
      r=2.0*M_PI*ran2(&seed);
#endif
      ps[i]=u*cos(r)+mu;
      ps[i+1]=u*sin(r)+mu;
   }     
   return 0;
}

References M_PI, ran2(), and seed.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ Gauss_f()

int Gauss_f	(	float *	ps,
		unsigned int	n,
		const float	mu,
		const float	sigma )

Generates a vector of normally distributed random single precision numbers using the Box-Muller Method.

Parameters

ps	The output array
n	The array length
mu	mean
sigma	variance

Returns: Zero on success integer error code otherwise

Definition at line 369 of file random.c.

                                                                         {
   /*
    * @brief   Generates a vector of normally distributed random single precision numbers using the Box-Muller Method
    * 
    * @param   ps:      The output array
    * @param   n:       The array length
    * @param   mu:      mean
    * @param   sigma:   variance
    *
    * @see ran2(), par_dcopy(), gsl_rng_uniform(), vdRngUniform()
    * 
    * @return Zero on success integer error code otherwise
    */
   const char *funcname = "Gauss_z";
   //The FORTRAN Code had two different Gauss Routines. gaussp having unit
   //mean and variance and gauss0 where the variance would appear to be 1/sqrt(2)
   //(Since we multiply by sqrt(-ln(r)) instead of sqrt(-2ln(r)) )
   if(n<=0){
      fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
            ARRAYLEN, funcname, n);
#if(nproc>1)
      MPI_Abort(comm,ARRAYLEN);
#else
      exit(ARRAYLEN);
#endif
   }
   int i;
   float r, u, v;
   //If n is odd we calculate the last index seperately and the rest in pairs
   if(n%2==1){
      n--;
#ifdef __RANLUX__
      r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
      ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
#else
      r=2.0*M_PI*ran2(&seed);
      ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
#endif
   }
#ifdef __RANLUX__
   r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
   ps[n]=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*cos(r);
#else
   r=2.0*M_PI*ran2(&seed);
   ps[n]=sqrt(-2*log(ran2(&seed)))*cos(r);
#endif
   for(i=0;i<n;i+=2){
      /* Marsaglia Method for fun
         do{
         u=sfmt_genrand_real1(sfmt);
         v=sfmt_genrand_real1(sfmt);
         r=u*u+v*v;
         }while(0<r & r<1);
         r=sqrt(r);
         r=sqrt(-2.0*log(r)/r)*sigma;
         ps[i] = mu+u*r; 
         ps[i+1]=mu+v*r;
       */
#ifdef __RANLUX__
      u=sqrt(-2*log(gsl_rng_uniform(ranlux_instd)))*sigma;
      r=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
#else
      u=sqrt(-2*log(ran2(&seed)))*sigma;
      r=2.0*M_PI*ran2(&seed);
#endif
      ps[i]=u*cos(r)+mu;
      ps[i+1]=u*sin(r)+mu;
   }
   return 0;
}

References M_PI, ran2(), and seed.

Here is the call graph for this function:

◆ Gauss_z()

int Gauss_z	(	Complex *	ps,
		unsigned int	n,
		const Complex	mu,
		const double	sigma )

Generates a vector of normally distributed random double precision complex numbers using the Box-Muller Method.

Parameters

ps	The output array
n	The array length
mu	mean
sigma	variance

Returns: Zero on success integer error code otherwise

Definition at line 214 of file random.c.

                                                                              {
   /*
    * @brief   Generates a vector of normally distributed random double precision complex numbers using the Box-Muller Method
    * 
    * @param   ps:      The output array
    * @param   n:       The array length
    * @param   mu:      mean
    * @param   sigma:   variance
    *
    * @see ran2(), par_dcopy(), gsl_rng_uniform(), vdRngUniform()
    * 
    * @return Zero on success integer error code otherwise
    */
   const char *funcname = "Gauss_z";
   if(n<=0){
      fprintf(stderr, "Error %i in %s: Array cannot have length %i.\nExiting...\n\n",
            ARRAYLEN, funcname, n);
#if(nproc>1)
      MPI_Abort(comm,ARRAYLEN);
#else
      exit(ARRAYLEN);
#endif
   }
#pragma unroll
   for(int i=0;i<n;i++){
      /* Marsaglia Method for fun
         do{
         u=sfmt_genrand_real1(sfmt);
         v=sfmt_genrand_real1(sfmt);
         r=u*u+v*v;
         }while(0<r & r<1);
         r=sqrt(r);
         r=sqrt(-2.0*log(r)/r)*sigma;
         ps[i] = mu+u*r + I*(mu+v*r);
       */
#ifdef __RANLUX__
      double   r =sigma*sqrt(-2*log(gsl_rng_uniform(ranlux_instd)));
      double   theta=2.0*M_PI*gsl_rng_uniform(ranlux_instd);
#else
      double   r =sigma*sqrt(-2*log(ran2(&seed)));
      double   theta=2.0*M_PI*ran2(&seed);
#endif
      ps[i]=r*(cos(theta)+sin(theta)*I)+mu;
   }     
   return 0;
}

References M_PI, ran2(), and seed.

Here is the call graph for this function:

◆ Par_granf()

double Par_granf ( )

Generates a random double which is then sent to the other ranks.

Returns: the random number generated

Definition at line 190 of file random.c.

                  {
   /*
    * @brief Generates a random double which is then sent to the other ranks
    *
    * @see ran2(), par_dcopy(), gsl_rng_uniform(), vdRngUniform()
    *
    * @return the random number generated
    */
   const char *funcname = "Par_granf";
   double ran_val=0;
   if(!rank){
#if (defined USE_RAN2||(!defined __INTEL_MKL__&&!defined __RANLUX__))
      ran_val = ran2(&seed);
#elif defined __RANLUX__
      ran_val = gsl_rng_uniform(ranlux_instd);
#else
      vdRngUniform(VSL_RNG_METHOD_UNIFORM_STD_ACCURATE, stream, 1, &ran_val, 0,1);
#endif
   }
#if(nproc>1)
   Par_dcopy(&ran_val);
#endif
   return ran_val;
}

References Par_dcopy(), ran2(), rank, and seed.

Here is the call graph for this function:

◆ Par_ranread()

int Par_ranread	(	char *	filename,
		double *	ranval )

Reads ps from a file Since this function is very similar to Par_sread, I'm not really going to comment it check there if you are confused about things.

Parameters

filename	The name of the file we're reading from
ranval	The destination for the file's contents

Returns: Zero on success, integer error code otherwise

Definition at line 88 of file random.c.

                                               {
   /*
    * @brief Reads ps from a file
    * Since this function is very similar to Par_sread, I'm not really going to comment it
    * check there if you are confused about things. 
    *
    * @param   *filename: The name of the file we're reading from
    * @param   ps:   The destination for the file's contents
    *
    * @return Zero on success, integer error code otherwise
    */
   const char *funcname = "Par_psread";
   FILE *dest;
   if(!rank){
      if(!(dest = fopen(filename, "rb"))){
         fprintf(stderr, "Error %i in %s: Failed to open %s.\nExiting...\n\n", OPENERROR, funcname, filename);
#if(nproc>1)
         MPI_Abort(comm,OPENERROR); 
#else
         exit(OPENERROR);
#endif
 
      }
      fread(&ranval, sizeof(ranval), 1, dest);  
      fclose(dest);
   }
#if(nproc>1)
   Par_dcopy(ranval);
#endif
   return 0;
}

References Par_dcopy(), and rank.

Here is the call graph for this function:

◆ Par_ranset()

int Par_ranset	(	long *	seed,
		int	iread )

Uses the rank to get a new seed. Copying from the FORTRAN description here c create new seeds in range seed to 9*seed c having a range of 0*seed gave an unfortunate pattern c in the underlying value of ds(1) (it was always 10 times bigger c on the last processor). This does not appear to happen with 9.

Parameters

seed	The seed from the rank in question.
iread	Do we read from file or not. Don't remember why it's here as it's not used

Returns: Zero on success, integer error code otherwise

Definition at line 135 of file random.c.

{
   const char *funcname = "Par_ranset";
   //If we're not using the master thread, we need to change the seed
#ifdef _DEBUG
   printf("Master seed: %i\t",*seed);
#endif
   if(rank)
      *seed *= 1.0f+8.0f*(float)rank/(float)(size-1);
#ifdef _DEBUG
   printf("Rank:  %i\tSeed %i\n",rank, *seed);
#endif
   //Next we set the seed using ranset
   //This is one of the really weird FORTRAN 66-esque functions with ENTRY points, so good luck!
#if (defined __INTEL_MKL__||defined __RANLUX__)
   return ranset(seed);
#else
   return 0;
#endif
}

References rank, ranset(), seed, and size.

Here is the call graph for this function:

Here is the caller graph for this function:

◆ ran2()

double ran2 ( long * idum )

Generates uniformly distributed random double between zero and one as described in numerical recipes. It's also thread-safe for different seeds.

Parameters

idum	Pointer to the seed

Returns: The random double between zero and one

Definition at line 440 of file random.c.

                        {
   /*
    * @brief   Generates uniformly distributed random double between zero and one as
    *          described in numerical recipes. It's also thread-safe for different seeds.
    *
    * @param   idum: Pointer to the seed
    *
    * @return  The random double between zero and one
    *
    */
   long k;
   int j;
   static long idum2=123456789; 
   static long iy=0;
   static long iv[NTAB];
#pragma omp threadprivate(idum2, iy, iv)
   //No worries
   double temp;
 
   if (*idum <= 0) {
      if (-(*idum) < 1) *idum=1; 
      else *idum = -(*idum);
      {
         idum2=(*idum);
 
         for(j=NTAB+7;j>=0;j--) {
            k=(*idum)/IQ1; 
            *idum=IA1*(*idum-k*IQ1)-k*IR1; 
            if (*idum < 0) *idum += IM1; 
            if (j < NTAB)
               iv[j] = *idum;
         }
         iy=iv[0];
      }
   }
   k=(*idum)/IQ1; 
   *idum=IA1*(*idum-k*IQ1)-k*IR1;
   if (*idum < 0) *idum += IM1; 
   k=idum2/IQ2; 
   idum2=IA2*(idum2-k*IQ2)-k*IR2;
 
   if (idum2 < 0) idum2 += IM2; j=iy/NDIV;
   iy=iv[j]-idum2;
   iv[j] = *idum;
   if (iy < 1) iy += IMM1;
   if ((temp=AM*iy) > RNMX) 
      return RNMX; 
 
   else return temp;
 
}

Here is the caller graph for this function:

◆ ranset()

int ranset ( long * seed )

inline

Dummy seed the ran2 generator.

Parameters

seed	pointer to seed

Returns: 0

Definition at line 74 of file random.c.

{
#ifdef __RANLUX__
   ranlux_instd=gsl_rng_alloc(gsl_rng_ranlxd2);
   gsl_rng_set(ranlux_instd,*seed);
   return 0;
#elif (defined __INTEL_MKL__&& !defined USE_RAN2)
   vslNewStream( &stream, VSL_BRNG_MT19937, *seed );
   return 0;
#else
   return 0;
#endif
}

References seed.

Here is the caller graph for this function:

Variable Documentation

◆ seed

long seed

extern

RAN2 seed.

Definition at line 27 of file random.c.

Macros

Functions

Variables

Detailed Description

Macro Definition Documentation

◆ AM

◆ EPS

◆ IA1

◆ IA2

◆ IM1

◆ IM2

◆ IMM1

◆ IQ1

◆ IQ2

◆ IR1

◆ IR2

◆ NDIV

◆ NTAB

◆ RAN2

◆ RNMX

Function Documentation

◆ Gauss_c()

◆ Gauss_d()

◆ Gauss_f()

◆ Gauss_z()

◆ Par_granf()

◆ Par_ranread()

◆ Par_ranset()

◆ ran2()

◆ ranset()

Variable Documentation

◆ seed