Monte Carlo estimation of π

Question

My C program uses the Monte Carlo method to approximate the mathematical constant π, the ratio of a circle's circumference to its diameter (and, importantly for this code, 4 times the ratio of a circle's area to that of its bounding square).

This estimation executes more slowly when parallelised using OpenMP than it does when executed in a single thread.

I think that my random generator function is not thread safe, and that may be the cause of the slowdown.

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

#define DART 1000000    // Number of darts each player throws
                        // Number of random dots in the square[-1,+1]

#define MAXPLAYER 8     // Maximume number of participant players
                    // Number of threads
/**
Generates a random number between two params given

@param: random number scope
@return: desired random number
*/
long double fRand(long double fMin, long double fMax)
{
    long double f = (double)rand() / RAND_MAX;
    return fMin + f * (fMax - fMin);
}

/**
    Generates random dots on the board

    @param: playerDarts, total number of darts to throw
    @return: score, number of darts thrown in the circle
*/ 
int player(int playersDarts)
{
    srand(time(NULL));
    long double pi, x, y;
    int score = 0;

    for (int i = 0; i < playersDarts; i++)
    {
        x = fRand(-1.0, 1.0);
        y = fRand(-1.0, 1.0);

        if (x*x + y*y < 1.0)
            score++;
    }
    return score;
}

void main()
{
long double pi;
long const double REAL_PI = 3.141592653589;
int score = 0, playersDarts;

////////////////////////////////////////////////////
//  Parallel                                      //
////////////////////////////////////////////////////

// devide the total number of DARTS between players
playersDarts = DART / MAXPLAYER;
double beginParallel = omp_get_wtime();

#pragma omp parallel for
for (int i = 1; i <= MAXPLAYER; i++)
    score += player(playersDarts);

double endParallel = omp_get_wtime();
pi = 4.0 * ((long double)score / (long double)DART);

double time_spent_parallel = endParallel - beginParallel;
printf("\n\t Calculated pi : %.12Lf\n", pi);
printf(  "\t       Real pi : %.12Lf\n", REAL_PI);
printf("\n\t Parallel Execution Time: %f\n", time_spent_parallel);

////////////////////////////////////////////////////
//  Serial                                        //
////////////////////////////////////////////////////

double beginSerial = omp_get_wtime();
pi = 4.0 * ((long double) player(DART) / (long double)DART);
double endSerial = omp_get_wtime();

double time_spent_serial = endSerial - beginSerial;

printf("\n\t Calculated pi : %.12Lf\n", pi);
printf(  "\t       Real pi : %.12Lf\n", REAL_PI);
printf("\n\t Serial Execution Time: %f\n", time_spent_serial);
}

Answer 1

In the code below you must realise that score is a variable that is shared between the threads. Which means that it requires synchronization, as you have omitted this synchornisation you will have undefined behaviour (incorrect result most likely).

for (int i = 1; i <= MAXPLAYER; i++)
    score += player(playersDarts);

OpenMP is not a magic box that will automatically make everything faster. You still need to think about synchronisation and proper algorithms. Also starting a thread does have some overhead and if the work done by the thread is small in comparison to the overhead, then you're slowing the program down.

There are plenty of resources on parallel programming on the internet, I recommend starting with doing your threading manually until you know how it works and then start using tools like OpenMP.

Answer 2

rand() is not thread safe. Calling it from multiple threads results in an undefined behavior. It may cause a slowdown, it may cause non-randomness, it may cause crash. You need rand_r() instead. See this answer for details.

Also notice that you never seed the RNG.

Answer 3

So according to your helpful answers and a little surfing on the web I found out that I need a pseudo random generator, which is a kind of RNG that does not generate repetitive numbers.

The other point is that I should have a different seed for each and every thread, so in the code I've passed the thread ID to the random generator function and used the rand_r() with different seeds.

here is the complete and no error code if you want to check it out: https://github.com/babak-ss/PiCalculationParallel

Answer 4

Simple mistakes:

• main() must return an int, and it's good practice to give it a full prototype: int main(void).
• Unused variable pi in player() function.

---

It's a bad idea to seed the random generator repeatedly from the current time (which changes very slowly compared to your program execution). Seed it once, around the beginning of the program (or once per thread; see below).

As others have noted, rand() is not thread-safe. We could create a critical section where you use it:

#pragma omp critical
    {
        x = fRand(-1.0, 1.0);
        y = fRand(-1.0, 1.0);
    }

However, that will really hammer performance here (where rand() is a large part of the runtime). We will be much better off using rand_r() (if you have it; it's required by POSIX but not standard C).

To use rand_r, we need a seed per thread, which needs to be passed through like this:

long double fRand(long double fMin, long double fMax, unsigned int *seedp)
{
    return fMin + (long double)rand_r(seedp) * (fMax - fMin) / RAND_MAX;
}


int player(int playersDarts)
{
    unsigned int seed;
#pragma omp critical
    {
        seed = (unsigned int)rand();
    }

    int score = 0;    
    for (int i = 0; i < playersDarts; i++)
        {
            long double x = fRand(-1.0, 1.0, &seed);
            long double y = fRand(-1.0, 1.0, &seed);

            if (x*x + y*y < 1.0)
                score++;
        }
    return score;
}

---

The += operation on score is not atomic, so it needs a critical section, or we can use an OpenMP reduction clause to give each thread its own copy, and add them all together at the end:

    int score = 0;
#pragma omp for reduction(+:score)
    for (int i = 0;  i < playersDarts;  ++i)
        {
            const long double x = fRand(-1.0, 1.0, &seed);
            const long double y = fRand(-1.0, 1.0, &seed);
            score += x*x + y*y < 1.0;
        }

---

Instead of having to specify the number of threads to use, we can just ask OpenMP to divide the work across all available cores (which won't suffer the possibility of rounding when DART isn't an exact multiple of MAXPLAYERS).

Here's a modified version that lets OpenMP divide the work amongst the cores - note how we can make the execution single-threaded just by setting OMP_NUM_THREADS in the environment:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <unistd.h>

#define DART 100000000    // Number of darts each player throws

long double fRand(long double fMin, long double fMax, unsigned int *seedp)
{
    return fMin + (long double)rand_r(seedp) * (fMax - fMin) / RAND_MAX;
}


int main(int argc, char **argv)
{
    long const double REAL_PI = 4 * atan(1);

    double beginTime, endTime;
    int score = 0;

    double beginTime = omp_get_wtime();
#pragma omp parallel
    {
        unsigned int seed;
#pragma omp critical
        {
            seed = (unsigned int)rand();
        }

#pragma omp for reduction(+:score)
        for (int i = 0;  i < DART;  ++i)
            {
                const long double x = fRand(-1.0, 1.0, &seed);
                const long double y = fRand(-1.0, 1.0, &seed);
                score += x*x + y*y < 1.0;
            }
    }
    double endTime = omp_get_wtime();

    long double pi = 4.0L * score / DART;

    printf("\n\t Calculated pi : %.12Lf\n", pi);
    printf(  "\t       Real pi : %.12Lf\n", REAL_PI);
    printf(  "\t         Error : %.3Lf%%\n", 100 * (pi - REAL_PI) / REAL_PI);
    printf("\n\t %d threads execution Time: %f\n", omp_get_max_threads(), endTime - beginTime);

    if (getenv("OMP_NUM_THREADS")) {
        return 0;
    }
    /* Do it all again, constrained to a single thread */
    setenv("OMP_NUM_THREADS", "1", 1);
    execvp(argv[0], argv);
}

I did have to add a couple of orders of magnitude to DARTS to get a meaningful difference in execution time between the single and parallel cases; the overhead of thread creation is simply too large otherwise, when the total runtime is under 100 ms.

I didn't seed the random number generator, so we get reproducible output. But if you do choose to do so, that must be before the #pragma omp parallel.