Monte Carlo pi逼近的并行化

我正在编写ac脚本来与OpenMp并行化pi近似。 我认为我的代码在令人信服的输出下运行良好。 我现在用4个线程运行它。 我不确定的是，如果此代码容易受到竞争条件的影响？ 如果是，如何协调此代码中的线程操作？

代码如下：

#include  #include  #include  #include  #include  double sample_interval(double a, double b) { double x = ((double) rand())/((double) RAND_MAX); return (ba)*x + a; } int main (int argc, char **argv) { int N = atoi( argv[1] ); // convert command-line input to N = number of points int i; int NumThreads = 4; const double pi = 3.141592653589793; double x, y, z; double counter = 0; #pragma omp parallel firstprivate(x, y, z, i) reduction(+:counter) num_threads(NumThreads) { srand(time(NULL)); for (int i=0; i < N; ++i) { x = sample_interval(-1.,1.); y = sample_interval(-1.,1.); z = ((x*x)+(y*y)); if (z<= 1) { counter++; } } } double approx_pi = 4.0 * counter/ (double)N; printf("%i %1.6e %1.6e\n ", N, 4.0 * counter/ (double)N, fabs(4.0 * counter/ (double)N - pi) / pi); return 0; } 

此外，我想知道是否应该在内部或外部并行化声明随机数的种子。 我的输出看起来像这样：

 10 3.600000e+00 1.459156e-01 100 3.160000e+00 5.859240e-03 1000 3.108000e+00 1.069287e-02 10000 3.142400e+00 2.569863e-04 100000 3.144120e+00 8.044793e-04 1000000 3.142628e+00 3.295610e-04 10000000 3.141379e+00 6.794439e-05 100000000 3.141467e+00 3.994585e-05 1000000000 3.141686e+00 2.971945e-05 

现在看起来还不错。 你对种族条件和种子安置的建议是最受欢迎的。

Answer 1

我可以看到您的代码中存在一些问题。 从我的观点来看，主要的是它没有并行化。 或者更准确地说，您在编译时没有启用OpenMP引入的并行性。 以下是人们可以看到的方式：

代码并行化的方式，主循环应该由所有线程完全执行（这里没有工作共享，没有#pragma omp parallel for ，只有#pragma omp parallel ）。 因此，考虑到将线程数设置为4，全局迭代次数应为4*N 因此，您的输出应该慢慢收敛到4 * Pi，而不是Pi。

实际上，我在我的笔记本电脑上尝试了你的代码，用OpenMP支持编译它，这就是我得到的。 但是，当我不启用OpenMP时，我会得到类似于你的输出。 总而言之，您需要：

1. 在编译时启用OpenMP以获取代码的并行版本。
2. 将您的结果除以NumThreads以获得Pi的“有效”近似值（或者例如使用#pragma omp for将您的循环分布在N上）

但是，如果你的代码在其他地方是正确的，那么它还没有。 正如BitTickler已经暗示的那样， rand()不是线程安全的。 所以你必须使用另一个随机数生成器，它允许你私有化它的状态。 例如，这可能是rand_r() 。 也就是说，这仍然有很多问题：

1. rand() / rand_r()在随机性和周期性方面是一个可怕的 RNG。 在增加尝试次数的同时，您将快速浏览RNG的周期并重复相同的序列。 你需要更强大的东西来做任何远程严肃的事情。
2. 即使使用“良好”的RNG，并行性方面也可能是一个问题，因为您希望并行的序列在彼此之间不相关。 并且每个线程使用不同的种子值并不能保证（虽然使用足够宽的RNG，你可以有一些余量）

无论如何，底线是：

• 使用更好的线程安全RNG（我发现drand48_r()或者random_r()对于Linux上的玩具代码没问题）
• 例如，根据线程id初始化其每个线程的状态，同时请记住，在某些情况下，这不会确保随机序列的正确去相关（并且调用函数的次数越多，可能性越大）你最终会有重叠的系列）。

这样做（以及一些小修复），您的代码变为例如如下：

 #include  #include  #include  #include  #include  typedef struct drand48_data RNGstate; double sample_interval(double a, double b, RNGstate *state) { double x; drand48_r(state, &x); return (ba)*x + a; } int main (int argc, char **argv) { int N = atoi( argv[1] ); // convert command-line input to N = number of points int NumThreads = 4; const double pi = 3.141592653589793; double x, y, z; double counter = 0; time_t ctime = time(NULL); #pragma omp parallel private(x, y, z) reduction(+:counter) num_threads(NumThreads) { RNGstate state; srand48_r(ctime+omp_get_thread_num(), &state); for (int i=0; i < N; ++i) { x = sample_interval(-1, 1, &state); y = sample_interval(-1, 1, &state); z = ((x*x)+(y*y)); if (z<= 1) { counter++; } } } double approx_pi = 4.0 * counter / (NumThreads * N); printf("%i %1.6e %1.6e\n ", N, approx_pi, fabs(approx_pi - pi) / pi); return 0; } 

我编译的是这样的：

 gcc -std=gnu99 -fopenmp -O3 -Wall pi.c -o pi_omp