The code below is a direct translation from a youtube video on Estimating PI using OpenMP and Monte Carlo. Even with the same inputs I'm not getting here their output. In fact, it seems like around half the value is what I get.

int main() {
    int num; // number of iterations
    printf("Enter number of iterations you want the loop to run for:   ");
    scanf_s("%d", &num);

    double x, y, z, pi;
    long long int i;
    int count = 0;
    int num_thread;

    printf("Enter number of threads you want to run to parallelize the process:\t");
    scanf_s("%d", &num_thread);
    printf("\n");

#pragma omp parallel firstprivate(x,y,z,i) shared(count) num_threads(num_thread)
    {
        srand((int)time(NULL) ^ omp_get_thread_num());
        for (i = 0; i < num; i  ) {
            x = (double)rand() / (double)RAND_MAX;
            y = (double)rand() / (double)RAND_MAX;

            z = pow(((x * x)   (y * y)), .5);

            if (z <= 1) {
                count  ;
            }
        }
    } // END PRAGMA
    pi = ((double)count / (double)(num * num_thread)) * 4;
    printf("The value of pi obtained is %f\n", pi);

    return 0;
}

I've also used a similar algorithm straight from the Oak Ridge National Laboratory's website (https://www.olcf.ornl.gov/tutorials/monte-carlo-pi/):

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>
 
int main(int argc, char* argv[])
{
    int niter = 1000000;            //number of iterations per FOR loop
    double x,y;                     //x,y value for the random coordinate
    int i;                          //loop counter
        int count=0;                //Count holds all the number of how many good coordinates
    double z;                       //Used to check if x^2 y^2<=1
    double pi;                      //holds approx value of pi
        int numthreads = 16;
 
    #pragma omp parallel firstprivate(x, y, z, i) shared(count) num_threads(numthreads)
    {
        srandom((int)time(NULL) ^ omp_get_thread_num());    //Give random() a seed value
        for (i=0; i<niter;   i)              //main loop
        {
            x = (double)random()/RAND_MAX;      //gets a random x coordinate
            y = (double)random()/RAND_MAX;      //gets a random y coordinate
            z = sqrt((x*x) (y*y));          //Checks to see if number is inside unit circle
        if (z<=1)
        {
                  count;            //if it is, consider it a valid random point
        }
    }
        //print the value of each thread/rank
    }
        pi = ((double)count/(double)(niter*numthreads))*4.0;
        printf("Pi: %f\n", pi);
 
    return 0;
}

And I am have the exact problem, so I'm think it isn't the code but somehow my machine.

I am running in VS Studio 22, Windows 11 with 16 core i9-12900kf and 32 gb ram.

Edit: I forgot to mention I did alter the second algorithm to use srand() and rand() instead.

CodePudding user response：

There are many errors in the code:

1. As pointed out by @JeromeRichard and @JohnBollinger rand\srand\random are not threadsafe you should use a threadsafe solution.

2. There is a race condition at line count; (different threads read and write a shared variable). You should use reduction to avoid it.

3. The code assumes that you use numthreads threads, but OpenMP does not guarantee that you actually got all of the threads you requested. I think if you got PI/2 as a result, the problem should be the difference between the requested and obtained number of threads. If you use #pragma omp parallel for... before the loop, you do not need any assumptions about the number of threads (ie. in this case the equation to calculate PI does not contain the number of threads).

4. A minor comment is that you do not need to use the time-consuming pow function.

Putting it together your code should be something like this:

#pragma omp parallel for reduction( :count) num_threads(num_thread)    
for (long long int i = 0; i < num; i  ) {
    const double x = threadsafe_random_number_between_0_1();
    const double y = threadsafe_random_number_between_0_1();
    const double z = x * x   y * y;
    if (z <= 1) {
        count  ;
    }
}
double pi = ((double) count / (double) num ) * 4.0;

CodePudding user response：

One assumption but I may be wrong : you initialise random with time, so it may happen than different thread use the same time , which may result in same random number generated, and so the result will be really bad as you got multiple time the same values. This is a problem with the Monte-Carlo method where 2 identical points will make wrong result.