Need help parallelising the ran2 and main program causing segmentation fault openmp

Question

I tried to parallelize the loop but increasing the thread count leave me with 'Segmentation fault(core dumped)'

I have parallelized the loop in the main.c, the loop internally refers to the ran2.c file function, which somewhere leads me to the segmentation fault, I need help parallelizing the program.

main.c:

#include <stdlib.h>
#include <stdio.h>
int main(int argc, char **argv) {
    int niter, i, j;
    long seed;
    double count;
    double x,y,z,pi;
    extern float ran2();
    niter=10000;
    count=0;

    #pragma omp parallel firstprivate(x, y, z, i) shared(count) num_threads(4)
    for(i=1; i<=niter; i++) {
        seed=i;
        x=ran2(&seed);
        y=ran2(&seed);
        z=x*x+y*y;
        if(z<1) {
            count+=1;
        }
    }
    pi=count*4.0/niter;
    printf("The value of pi is %8.14f\n",pi);
    return 0;
}

ran2.c

#define IM2 2147483399
#define IM1 2147483563
#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)

float ran2(long *idum) {
    int j;
    long k;
    static long idum2=123456789;
    static long iy=0;
    static long iv[NTAB];
    float 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;
}

I need to parallelize the program.

Answer 1

Your code has several problems.

First of all, you are passing the uninitialized values of x, y, and z to each OpenMP thread with firstprivate, which the compiler warns about. This is easy to work around by moving the declarations inside to for loop and removing the firstprivate. For the same reason, i should be private.

seed is supposed to stay the same during each iteration, but since it is a shared variable, the threads are competing to get their value written to it. This is also fixed by moving the declaration of seed into the loop.

The third problem is that the threads are also racing when writing to count. Since you want to sum up the values from the different threads, you should be using the reduction(+:count) directive.

Fourth, the ran2 function is not thread-safe. It has static state in the form of idum2, iy, and iv, where again the threads race. The best way to fix this is to get rid of the global state by moving it into a struct ran_state (see below) which is passed as a parameter. Then replace each place referring to idum2 with state->idum2, iy with iy and state->iv with state->iv. Finally, create and initialize this state before the for loop, and mark it as firstprivate so each thread gets its own initialized state.

This then looks as below (in one file for convenience):

#include <stdlib.h>
#include <stdio.h>

// Insert the #defines for ran2 here ...

struct ran_state {
    long idum2;
    long iy;
    long iv[NTAB];
};

// Insert ran2 here, replacing iv with state->iv, and so on ...

int main(int argc, char **argv){
    int niter, i;
    double count;
    double pi;
    niter=10000;
    count=0;
    struct ran_state state = { 123456789, 0, {0} };

#pragma omp parallel for reduction(+:count) firstprivate(state) private(i) num_threads(4)
    for(i=1;i<=niter;i++){
        double x, y, z;
        long seed;
        seed=i;
        x=ran2(&seed, &state);
        y=ran2(&seed, &state);
        z=x*x+y*y;
        if(z<1){
            count+=1;
        }
    }
    pi=count*4.0/niter;
    printf("The value of pi is %8.14f\n",pi);
    return 0;
}