Random number generator related routines. More...

#include "coord.h"
#include "errorcodes.h"
#include "par_mpi.h"
#include "random.h"
#include "sizes.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>

Include dependency graph for random.c:

Macros
#define	M_PI acos(-1)
	\(\pi\) if not defined elsewhere

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

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.

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	Par_granf ()
	Generates a random double which is then sent to the other ranks.

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_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_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_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.

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.

Variables
long	seed
	RAN2 seed.

Detailed Description

Random number generator related routines.

Definition in file random.c.

Macro Definition Documentation

◆ M_PI

#define M_PI acos(-1)

\(\pi\) if not defined elsewhere

Definition at line 41 of file random.c.

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

RAN2 seed.

Definition at line 27 of file random.c.

Macros

Functions

Variables

Detailed Description

Macro Definition Documentation

◆ M_PI

Function Documentation

◆ Gauss_c()

◆ Gauss_d()

◆ Gauss_f()

◆ Gauss_z()

◆ Par_granf()

◆ Par_ranread()

◆ Par_ranset()

◆ ran2()

◆ ranset()

Variable Documentation

◆ seed