如何使用 pthread 并行化 pi 的蒙特卡罗估计？

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 );
}

上面的代码是一个顺序程序，它接受一个参数来表示估计 pi ​​的接近程度。一旦这些新旧值的变化低于它退出的容差。

我必须在 pthreads 中并行化这个程序。我不是想让某人为我做这件事，而是想得到一些建议和想法来思考，以便我能够做到这一点。 pthreads 程序将线程数和容差作为它的参数并输出估计值。我对并行程序很陌生，真的不知道从哪里开始，所以我会接受任何建议。谢谢。

Answer 1

每个线程都应该保持自己的 in 和总计数，并使用更好的退出条件。然后将每个线程的 in 和 total 值相加。

Answer 2

#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;
}

这是我突然想到的，但我遇到了错误。它只是停止。我只是让线程跳出循环并返回，然后主线程加入它们。对我在这里做什么有什么建议吗？

我运行它，似乎所有线程在达到互斥锁时都被锁定了。我让每个线程每 1000 点检查一次 pi 的变化。