How to parallelize monte carlo estimation of pi using pthreads?

Question

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

int main(int argc, char **argv)
{

    unsigned long long in = 1;
    unsigned long long total = 2;
    double tol , change, new, secs , old = 0.0;
    struct timeval start , end;
    int threads ; /* ignored */

    if ( argc < 2) {
        exit (-1);
    }

    threads = atoi ( argv[1]);
    tol = atof ( argv[2]);
    if (( threads < 1) || ( tol < 0.0)) {
        exit (-1);
    }
    tol = tol *tol;

    srand48(clock());
    gettimeofday (&start , NULL);
    do
    {
        double x, y;
        x = drand48();
        y = drand48();
        total ++;
        if (( x*x + y*y) <= 1.00)
            in ++;
        new = 4.0 * (double)in/( double)total ;
        change = fabs (new - old);
        old = new;
    }while (change > tol );
    gettimeofday (&end, NULL);
    secs = (( double)end.tv_sec - (double)start.tv_sec )
    + (( double)end.tv_usec - (double)start.tv_usec )/1000000.0;
    printf ( ”Found estimate of pi of %.12f in %llu iterations , %.6f seconds.n n”,
            new, total - 2, secs );
}

The code above is a sequential program that takes an argument for the tolerance for how closely to estimate pi. Once the change in these old and new values fall below the tolerance it exits.

I have to parallelize this program in pthreads. I am not trying to have someone do it for me but rather to get some pointers and ideas to think about so that I may be able to this. the pthreads program will take number of threads and tolerance as it's argument and output the estimation. I am very new to parallel programs and don't really know where to start so I will take any advice. Thanks.

Answer 1

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <unistd.h>
#include <pthread.h>

void* Function(void* i);

#define MAX_THREADS 200
unsigned long long total[MAX_THREADS] = {0};    //total points for threads
unsigned long long in[MAX_THREADS] = {0};       //points inside for threads
double tolerance, change, new, old = 0.0;
long thread_num;
pthread_mutex_t         mutex = PTHREAD_MUTEX_INITIALIZER;



int main(int argc, char **argv)
{
    long i;
    struct timeval start, end;
    double secs;
    unsigned long long master_total;
    pthread_t* threads;

    if (argc != 3){
        printf("\nMust pass 2 arguments:  (Tolerance) (# of Threads)");
        exit(-1);
    }

    thread_num = atoi ( argv[1]);
    tolerance = atof ( argv[2]);

    if (( thread_num < 1) || ( tolerance < 0.0) || (thread_num > 200)) {
        printf("\nIncorrect tolerance or threads.");
        exit (-1);
    }

    threads = malloc(thread_num*sizeof(pthread_t)); //allocating space for threads
    tolerance = tolerance * tolerance;
    change = 0.5;
    srand48(clock());
    gettimeofday (&start, NULL);
    for( i = 0; i < thread_num; i++ ){
        pthread_create(&threads[i], NULL, Function, (void*)i);
    }

    for( i = 0; i < thread_num; i++ ){
        pthread_join(threads[i], NULL);
    }
    gettimeofday (&end, NULL);

    master_total = 0;
    for( i = 0; i < thread_num; i++ ){
        master_total = master_total + total[i];
    }
    secs = (( double)end.tv_sec - (double)start.tv_sec )
    + (( double)end.tv_usec - (double)start.tv_usec )/1000000.0;
    printf ( "Estimation of pi is %.12f in %llu iterations , %.6f seconds.\n", new, master_total, secs );

}
//Each thread will calculate it's own points for in and total
//Per 1000 points it will calculate new and old values and compare to tolerance
//If desired accuracy is met the threads will return. 
void* Function(void* i){
    /*
     rc - return code
     total[i], in[i] - threads own number of calculated points
     my_total, my_in - Each thread calculates global in and total, per 1000 points and calculates tolerance

     */
    long my_spot = (long) i;
    long rc;
    long j;
    unsigned long long my_total;
    unsigned long long my_in;
    do
    {
        double x, y;
        x = drand48();
        y = drand48();
        total[my_spot]++;
        if (( x*x + y*y) <= 1.00){
            in[my_spot]++;
        }
        if(total[my_spot] % 1000 == 0){
            while ( j < thread_num){
                my_total = my_total + total[j];
                my_in = my_in + in[j];
            }
            my_total = my_total;
        //critical section
        //old, new, and change are all global
        rc = pthread_mutex_lock(&mutex);
        new = 4.0 * (double)my_in/( double)my_total;
        change = fabs (new - old);
        old = new;
        rc = pthread_mutex_unlock(&mutex);
        }
    }while (change > tolerance );
    return NULL;
}

This is what I whipped up but I am getting errors. It just stops. I just have the threads break out of a loop and return then the main thread joins them. Any advice on what I am doing here?

I run it and it seems that all the threads get locked out when they reach mutex lock. I have each thread check the change in pi every 1000 points.

Answer 2

Each thread should keep its own in and total counts and use a better exit condition. Then add up the in and total values from each thread.