用非线程安全随机数生成器修正C语言中pi monte carlo的OpenMP语用

Question

我需要一些帮助，用一个给定的随机数生成器用openmp并行化与蒙特卡洛方法的pi计算，这不是线程安全的。

第一：这条线对我没有帮助。

我自己的尝试是下面的#pragma omp语句。我认为I、x和y变量应该是每个线程的init，并且应该是private的。z是圆圈中所有命中数的总和，因此它应该在for循环后隐含的barriere之后求和。

认为主要的问题是随机数发生器的静态无功。我在调用函数的地方做了一个关键部分，这样每次只有一个线程可以执行它。但Pi解并不具有更高的值。

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

    int i, z = 0, threads = 8, iters = 100000;
    double x,y, pi;

    #pragma omp parallel firstprivate(i,x,y) reduction(+:z) num_threads(threads)
        for (i=0; i<iters; ++i) {
            #pragma omp critical
            {
                x = rng_doub(1.0);
                y = rng_doub(1.0);
            }
            if ((x*x+y*y) <= 1.0)
                z++;
        }

    pi = ((double) z / (double) (iters*threads))*4.0;
    printf("Pi: %lf\n", pi);;
    return 0;
}

#define RNG_MOD 741025

int rng_int(void) {
    static int state = 0;

    return (state = (1366 * state + 150889) % RNG_MOD);
}

double rng_doub(double range) {
    return ((double) rng_int()) / (double) ((RNG_MOD - 1)/range);
}

我也尝试过使静态int状态全局化，但它并没有改变我的结果，也许我做错了。所以你能帮我做正确的修改吗？非常感谢！

Answer 1

请尝试下面的代码。它为每个线程创建一个私有状态。我对atrand_r函数做了类似的操作，为什么使用OpenMP进行计算要比使用单个线程多花100倍的时间呢？

编辑：我使用Hristo的一些建议更新了我的代码。我（第一次）使用了threadprivate。我还使用了一个更好的兰德函数，它给出了pi的更好估计，但它仍然不够好。

一个奇怪的事情是我必须在threadprivate之后定义函数rng_int否则我会得到一个错误“error：'state'声明'threadprivate'after first use”。我也许该问个问题。

//gcc -O3 -Wall -pedantic -fopenmp main.c
#include <omp.h>
#include <stdio.h>

#define RNG_MOD 0x80000000
int state;

int rng_int(void);
double rng_doub(double range);

int main() {
    int i, numIn, n;
    double x, y, pi;

    n = 1<<30;
    numIn = 0;

    #pragma omp threadprivate(state)
    #pragma omp parallel private(x, y) reduction(+:numIn) 
    {

        state = 25234 + 17 * omp_get_thread_num();
        #pragma omp for
        for (i = 0; i <= n; i++) {
            x = (double)rng_doub(1.0);
            y = (double)rng_doub(1.0);
            if (x*x + y*y <= 1) numIn++;
        }
    }
    pi = 4.*numIn / n;
    printf("asdf pi %f\n", pi);
    return 0;
}

int rng_int(void) {
   // & 0x7fffffff is equivalent to modulo with RNG_MOD = 2^31
   return (state = (state * 1103515245 + 12345) & 0x7fffffff);
}

double rng_doub(double range) {
    return ((double)rng_int()) / (((double)RNG_MOD)/range);
}

您可以在http://coliru.stacked-crooked.com/a/23c1753a1b7d1b0d上看到结果（以及编辑和运行代码