Monte Carlo pi 近似的并行化

Question

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

代码如下所示：

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <omp.h>

double sample_interval(double a, double b) {

  double x = ((double) rand())/((double) RAND_MAX);
  return (b-a)*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 时没有启用您在 OpenMP 中引入的并行性。这是人们可以看到的方式：

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

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

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

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

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

无论如何，底线是：

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

完成此操作（以及一些小修复），您的代码例如如下所示：

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include <math.h>
#include <omp.h>

typedef struct drand48_data RNGstate;

double sample_interval(double a, double b, RNGstate *state) {
    double x;
    drand48_r(state, &x);
    return (b-a)*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