Monte Carlo pi 近似的并行化
Parallelization for Monte Carlo pi approximation
我正在编写一个 c 脚本来将 pi 近似与 OpenMp 并行化。我认为我的代码运行良好,输出令人信服。我现在 运行 它有 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
目前看起来还不错。非常欢迎您提出有关比赛条件和种子放置的建议。
我看到您的代码中存在一些问题。主要的是从我的角度来看,它不是并行化的。或者更准确地说,您在编译时没有启用您在 OpenMP 中引入的并行性。这是人们可以看到的方式:
代码并行化的方式,主for
循环应该被所有线程完整执行(这里没有工作共享,没有#pragma omp parallel for
,只有一个#pragma omp parallel
).因此,考虑到你设置的线程数为4,全局迭代次数应该是4*N
。因此,您的输出应该慢慢收敛到 4*Pi,而不是 Pi。
的确,我在我的笔记本电脑上试过你的代码,在 OpenMP 支持下编译它,这几乎就是我得到的。但是,当我不启用 OpenMP 时,我会得到类似于您的输出。所以总而言之,您需要:
- 在编译时启用 OpenMP 以获得代码的并行版本。
- 将你的结果除以
NumThreads
得到 Pi 的 "valid" 近似值(或者将你的循环分布在 N
上 #pragma omp for
例如)
但那是如果/当您的代码在其他地方是正确的,但现在还不是。
正如 BitTickler 已经暗示的那样,rand()
不是线程安全的。所以你必须去寻找另一个随机数生成器,这将允许你私有化它的状态。例如,这可能是 rand_r()
。也就是说,这仍然有很多问题:
rand()
/ rand_r()
在随机性和周期性方面是 糟糕的 RNG。在增加尝试次数的同时,您将快速遍历 RNG 期间并一遍又一遍地重复相同的序列。你需要更强大的东西来做任何严肃的事情。
- 即使使用 "good" 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
我正在编写一个 c 脚本来将 pi 近似与 OpenMp 并行化。我认为我的代码运行良好,输出令人信服。我现在 运行 它有 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
目前看起来还不错。非常欢迎您提出有关比赛条件和种子放置的建议。
我看到您的代码中存在一些问题。主要的是从我的角度来看,它不是并行化的。或者更准确地说,您在编译时没有启用您在 OpenMP 中引入的并行性。这是人们可以看到的方式:
代码并行化的方式,主for
循环应该被所有线程完整执行(这里没有工作共享,没有#pragma omp parallel for
,只有一个#pragma omp parallel
).因此,考虑到你设置的线程数为4,全局迭代次数应该是4*N
。因此,您的输出应该慢慢收敛到 4*Pi,而不是 Pi。
的确,我在我的笔记本电脑上试过你的代码,在 OpenMP 支持下编译它,这几乎就是我得到的。但是,当我不启用 OpenMP 时,我会得到类似于您的输出。所以总而言之,您需要:
- 在编译时启用 OpenMP 以获得代码的并行版本。
- 将你的结果除以
NumThreads
得到 Pi 的 "valid" 近似值(或者将你的循环分布在N
上#pragma omp for
例如)
但那是如果/当您的代码在其他地方是正确的,但现在还不是。
正如 BitTickler 已经暗示的那样,rand()
不是线程安全的。所以你必须去寻找另一个随机数生成器,这将允许你私有化它的状态。例如,这可能是 rand_r()
。也就是说,这仍然有很多问题:
rand()
/rand_r()
在随机性和周期性方面是 糟糕的 RNG。在增加尝试次数的同时,您将快速遍历 RNG 期间并一遍又一遍地重复相同的序列。你需要更强大的东西来做任何严肃的事情。- 即使使用 "good" 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