OpenMP Monte_Carlo симуляция достижения целевой близости к PI - PullRequest
Новая
       11

OpenMP Monte_Carlo симуляция достижения целевой близости к PI

0 голосов
/ 24 ноября 2018

Я пытаюсь написать параллельную программу, которая принимает коэффициент ошибок (т. Е. 0,01) и возвращает значение PI, которое ближе к PI, чем ошибка с симуляцией Монте-Карло.Я написал простую функцию, однако она не заканчивается, так как частота ошибок всегда составляет около 11. Я ценю ваши комментарии. 

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

double drand48(void);

double monte_carlo(double epsilon){
    double x,y, pi_estimate = 0.0;
    double drand48(void);
    double error = 10000.0;
    int n = 0; // total number of points
    int i = 0; // total numbers of points inside circle
    int p = omp_get_num_threads();
    while(error>=epsilon){
        #pragma omp parallel private(x, y) reduction(+:i)//OMP parallel directive
        {
            x = drand48();
            y = drand48();
            if((x*x+y*y)<=1.0){i+=1;}
        }
        n+=p;
        printf("%lf\n", error);
        pi_estimate=4.0*(double)i/(double)n;
        error = fabs(M_PI-pi_estimate)/M_PI;
    }
    return pi_estimate;
}

int main(int argc, char* argv[]) {
    double epsilon = 0.01;
    printf("PI estimate: %lf",monte_carlo(epsilon));
    return 0;
}

1 Ответ

0 голосов
/ 26 ноября 2018

Вызов omp_get_num_threads() вне параллельной секции всегда вернет 1, так как в данный момент вызывается только один активный поток.Следующий код должен дать правильный результат, но он будет намного медленнее, чем последовательная версия, из-за больших затрат на параллелизацию и синхронизацию для выполнения очень простой операции. 

#pragma omp parallel private(x, y) reduction(+:i)//OMP parallel directive
{
    x = drand48();
    y = drand48();
    if((x*x+y*y)<=1.0){i+=1;}
    #pragma omp master
    n+=omp_get_num_threads();
}

Следующее позволяет избежать многократного порождения потоков иможет быть более эффективным, но все же, вероятно, медленнее. 

#pragma omp parallel private(x, y)
while(error>=epsilon){
        x = drand48();
        y = drand48();
        if((x*x+y*y)<=1.0){
            #pragma omp atomic
            i++;
        }
    #pragma omp barrier
    #pragma omp single
    {
        n+=omp_get_num_threads();
        pi_estimate=4.0*(double)i/(double)n;
        error = fabs(M_PI-pi_estimate)/M_PI;
        printf("%lf\n", error);
    } // implicit barrier here
}

Для того, чтобы действительно идти быстрее, нужно указать минимальное количество итераций, например: 

#define ITER 1000
#pragma omp parallel private(x, y)
while(error>=epsilon){
    #pragma omp for reduction(+:i)
    for (int j=1;j<ITER;j++){
        x = drand48();
        y = drand48();
        if((x*x+y*y)<=1.0) i+=1;
    }
    /* implicit barrier + implicit atomic addition
     * of thread-private accumulator to shared variable i
     */

    #pragma omp single
    {
        n+=ITER;
        pi_estimate=4.0*(double)i/(double)n;
        error = fabs(M_PI-pi_estimate)/M_PI;
        printf("%lf\n", error);
    } // implicit barrier
}
...