c - 并行 block 中允许哪些功能和操作？

Question

代码:

double x(){return (double)rand()/(double)RAND_MAX;}
double y(){return (double)rand()/(double)RAND_MAX;}
double z(){return (double)rand()/(double)RAND_MAX;}

int d(double x, double y, double z){
        if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
        return 0;
    }

double f(double x, double y, double z){
        return 1;
    }




#pragma omp parallel default(none) private(id,numt,j,local_sum,local_good_dots,local_coi,x_,y_,z_) shared(total_sum,good_dots,count_of_iterations)
    {
        local_coi = count_of_iterations;
        id = omp_get_thread_num() + 1;
        numt = omp_get_num_threads();
        #pragma omp for
        for (j = 1; j <= local_coi;  j++){
            x_=x();
            y_=y();
            z_=z();
            if (d(x_,y_,z_) == 1){
                local_sum += f(x_,y_,z_);
                local_good_dots += 1;

            }
        }

        #pragma omp critical
        {
            total_sum = total_sum + local_sum;
            good_dots = good_dots + local_good_dots;
        }
    }

注释:该代码是蒙特卡罗方法的实现，用于计算函数f()在区域d()中的三维积分。

我希望这段代码在多线程模式(openmp)下运行得更快。

但是出了点问题。

经过几个小时的修改(openmp pragma 中的减少、if 条件的简化(例如 f(x_,y_,z_) * d(x_,y_,z_) code>)) 我不明白，为什么这个简单的循环在更多线程上变得更慢。

但是当我在循环之前为每个坐标生成一个 3 维数组并将其放入 shared 后，我的程序变得更快。

那么，问题:

如何修改此代码以及并行 block 中允许哪些功能(操作)？

P.S:据我所知，不允许使用 rand 函数(或者我错了？)

感谢您的帮助!

修改(在 @HristoIliev 的帮助下)

double x(){return (double)rand()/(double)RAND_MAX;}
double y(){return (double)rand()/(double)RAND_MAX;}
double z(){return (double)rand()/(double)RAND_MAX;}

int d(double x, double y, double z){
        if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
        return 0;
    }

double f(double x, double y, double z){
        return 1;
    }


#pragma omp parallel default(none) private(j,local_coi,x_,y_,z_) shared(count_of_iterations) reduction(+:total_sum,good_dots)
    {
        local_coi = count_of_iterations;
        #pragma omp for(prng)
        for (j = 1; j <= local_coi;  j++){                    
        #pragma omp critical(prng)
        {
                x_=x();
                y_=y();
                z_=z();
        }   
            if (d(x_,y_,z_) == 1){
                total_sum += f(x_,y_,z_);
                good_dots += 1;

            }
        }
    }

Answer 1

随机数生成器rand()使用全局静态分配状态，由所有线程共享，因此不是线程安全的。从多个线程使用它，您会遇到对共享变量进行不 protected 访问的非常糟糕的情况，这会破坏缓存并减慢程序速度。您应该使用 rand_r() 或 erand48() - 它们使用您必须提供的单独状态存储。您必须为每个线程声明一个状态(例如将其设为私有(private))，基本上为每个线程创建不同的 PRNG。然后你必须相应地播种它们，否则你会得到统计上不好的结果。原则上，您可以使用一个 rand48() 生成器的输出来为其他生成器提供种子 - 它应该足以获得中等长度的不相关序列。

这里是使用 rand_r() 的示例实现(并不是说这是一个用于蒙特卡洛模拟的非常糟糕生成器，erand48 更好最好是使用 GNU Scientific Library 中的“Mersenne Twister”类型生成器(如果可用):

unsigned int prng_state;
#pragma omp threadprivate(prng_state)

double x(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double y(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}
double z(){return (double)rand_r(&prng_state)/(double)RAND_MAX;}

int d(double x, double y, double z){
    if ( ( (pow(x,2)+pow(y,2)) <1 ) && ( z<=1 && z>=0 )) return 1;
    return 0;
}

double f(double x, double y, double z){
    return 1;
}

...

#pragma omp parallel default(none) \
            private(id,numt,x_,y_,z_) \
            shared(count_of_iterations) \
            reduction(+:total_sum,good_dots)
{
    id = omp_get_thread_num() + 1;
    numt = omp_get_num_threads();

    // Sample PRNG seeding code - DO NOT USE IN PRODUCTION CODE!
    prng_state = 67894 + 1337*id;

    #pragma omp for
    for (j = 1; j <= count_of_iterations;  j++){
        x_=x();
        y_=y();
        z_=z();
        if (d(x_,y_,z_) == 1){
            total_sum += f(x_,y_,z_);
            good_dots += 1;
        }
    }
}

这只是一个非常糟糕的(从质量角度来看)实现，但它应该让您了解事情是如何工作的。这也是您如何通过对原始代码进行最少的更改来实现线程安全。基本要点是:

• OpenMP threadprivate 指令将 PRNG 状态 prng_state 设置为每个线程私有(private)；
在 x()、 y 中，使用具有线程特定状态变量的
• rand_r() 代替 rand() () 和 z();
• PRNG 状态以线程相关的方式初始化，例如prng_state = 67894 + 1337*id;，以便不同的线程(希望)获得不相关的伪随机数流。

请注意，rand() 和 rand_r() 的质量很差，这只是一个学术示例。使用较长的 PRNG 序列，您将在不同线程中获得相关流，这会破坏统计数据。我让您使用 erand48() 重写代码。

回答您最初的问题 - 所有线程安全函数调用都允许在并行 block 内。您还可以调用非线程安全函数，但必须保护(命名)关键构造内的调用，例如:

#pragma omp for
for (j = 1; j <= local_coi; j++) {
    #pragma omp critical(prng)
    {
        x_=x();
        y_=y();
        z_=z();
    }
    if (d(x_,y_,z_) == 1) {
        local_sum += f(x_,y_,z_);
        local_good_dots += 1;
    }
}

这将确保不会并行调用rand()。但您仍然可以对共享状态进行读-修改-写访问，因此会出现与缓存相关的速度减慢。

此外，请勿尝试重新实现 OpenMP reduction 或类似构造。编译器供应商已经投入巨大的努力来确保它们以尽可能最好(读取最快)的方式实现。