OpenMP C Normal Random Generator

Question

I was doing a C assignment for parallel computing, where I have to implement some sort of Monte Carlo simulations with efficient tread safe normal random generator using Box-Muller transform. I generate 2 vectors of uniform random numbers X and Y, with condition that X in (0,1] and Y in [0,1]. But I'm not sure that my way of sampling uniform random numbers from the halfopen interval (0,1] is right. Did anyone encounter something similar? I'm using following Code:

double* StandardNormalRandom(long int N){
    double *X = NULL, *Y = NULL, *U = NULL;

    X = vUniformRandom_0(N / 2);
    Y = vUniformRandom(N / 2);
    #pragma omp parallel for
    for (i = 0; i<N/2; i++){
            U[2*i] = sqrt(-2 * log(X[i]))*sin(Y[i] * 2 * pi);
            U[2*i + 1] = sqrt(-2 * log(X[i]))*cos(Y[i] * 2 * pi);
            }
    return U;
    }

    double* NormalRandom(long int N, double mu, double sigma2)
    {
            double *U = NULL, stdev = sqrt(sigma2);
            U = StandardNormalRandom(N);
            #pragma omp parallel for
            for (int i = 0; i < N; i++) U[i] = mu + stdev*U[i];
            return U;
    }

here is the bit of my UniformRandom function also implemented in parallel:

#pragma omp parallel for  firstprivate(i)
for (long int j = 0; j < N;j++)
{
    if (i == 0){
        int tn = omp_get_thread_num();
        I[tn] = S[tn];
        i++;
    }
    else
    {
        I[j] = (a*I[j - 1] + c) % m;
    }

}

}
#pragma omp parallel for 
for (long int j = 0; j < N; j++)
    U[j] = (double)I[j] / (m+1.0);

score 0 · Accepted Answer · answered Jan 10 '14 at 20:14

0

In the StandardNormalRandom function, I will assume that the pointer U has been allocated to the size N, in which case this function looks fine to me. As well as the function NormalRandom.

However for the function UniformRandom (which is missing some parts, so I'll have to assume some stuff), if the following line I[j] = (a*I[j - 1] + c) % m + 1; is the body of a loop with a omp parallel for, then you will have some issues. As you can't know the order of execution of the thread, the current thread (with a fixed value of j) can't rely on the value of I[j - 1] as this value could be modified at any time (I should be shared by default).

Hope it helps!

answered Jan 10 '14 at 20:14

A. Lefebvre

139
4

This is not an issue, I initialise I[j] at the first execution of the thread and give it value from the array of random numbers I generated before. – Suren Jan 12 '14 at 13:02
@Suren Ok, could you edit your question with the missing code. It will be easier for people to help you. – A. Lefebvre Jan 13 '14 at 14:31

OpenMP C Normal Random Generator

1 Answers1