0

代码：

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）下运行得更快。

但是出了点问题。

经过几个小时的修改（reduction在 openmp pragma 中，if-condition 的简化（如f(x_,y_,z_) * d(x_,y_,z_)））我不明白，为什么这个简单的循环在更多的线程上变得更慢。

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

所以，问题：

如何修改此代码以及并行块中允许哪些功能（操作）？

PS：如我所见，该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;

            }
        }
    }
4

1 回答 1

2

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

这是一个使用示例实现rand_r()（并不是说这是一个非常糟糕的蒙特卡洛模拟生成器，erand48更好，最好是使用来自 GNU 科学库的“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指令，PRNG 状态prng_state对每个线程都是私有的；threadprivate
  • rand_r()使用线程特定的状态变量代替rand()in x(),y()和z();
  • PRNG 状态以与线程相关的方式初始化，例如prng_state = 67894 + 1337*id;，以便不同的线程（希望）获得不相关的伪随机数流。

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

要回答您最初的问题 - 块内允许所有线程安全的函数调用parallel。您也可以调用非线程安全函数，但必须保护（命名）critical构造内部的调用，例如：

#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()并行进行调用。但是您仍然可以对共享状态进行读-修改-写访问，因此与缓存相关的减速。

此外，不要尝试重新实现 OpenMPreduction或类似结构。编译器供应商已经付出了巨大的努力来确保它们以尽可能最好（读取速度最快）的方式实现。

于 2012-11-19T17:08:14.860 回答