重要 TL;DR rand is not thread safe:
来自rand 手册页:
函数 rand() 不是可重入的或线程安全的,因为它使用在每次调用时都会修改的隐藏状态。
对于多线程代码使用(例如)rand_r。
我试图了解 OpenMP 中的缩减是如何工作的。
为了论证起见,我们假设r2() 将始终产生相同的值。
当一个人的代码有多个线程同时修改某个变量时,代码如下所示:
double S = 0;
#pragma omp parallel
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
变量S 的更新有一个竞争条件。为了解决这个问题,可以使用 OpenMP reduction 子句,从 OpenMP standard 可以读取:
归约子句可用于执行某些形式的递归
并行计算 (...)。 用于并行和工作共享
构造,创建每个列表项的私有副本,每个列表项一个
隐式任务,好像使用了 private 子句。 (...) 这
然后按上面指定的方式初始化私有副本。在结束时
指定缩减子句的区域,原始列表
通过将其原始值与最终值组合来更新项目
每个私有副本,使用指定的组合器
减少标识符。
在这种情况下,代码如下所示:
#pragma omp for reduction(+:S)
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
但是,在您的完整代码中
#pragma omp parallel for collapse(2)
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
#pragma omp for reduction(+:S)
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
K[m * N + n] = S;
}
}
您首先使用#pragma omp for collapse(2) 划分两个外部循环的迭代,然后尝试使用不同的子句#pragma omp for 再次划分最内部循环的迭代,这是不允许的。
这是执行嵌套循环的正确方法吗?
您可以进行以下并行化:
#pragma omp parallel for collapse(2) firstprivate (S)
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
K[m * N + n] = S;
}
}
没有竞争条件,因为变量S 是私有的。此外,在这种情况下,由于两个最外层循环的迭代在线程之间划分,每个线程都有一对唯一的m和n迭代,因此每个线程将访问数组K的唯一位置。它的访问权限K[m * N + n]。
但问题是并行化两个外部循环的版本不会产生与其顺序对应的相同结果。之所以如此,是因为
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
K[m * N + n] = S;
在三个循环的所有迭代中添加一个隐式依赖。
S 的值明确取决于迭代 m、n 和 o 的执行顺序。因此,如果在线程之间划分这些循环的迭代,那么如果代码按顺序或并行执行,则给定m 和n 的S 的值将不同。尽管如此,这可以通过仅并行化最内层循环并减少变量S来解决:
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
#pragma omp parallel for reduction(+:S)
for (int o = 0; o < O; o++) {
S += r2() - 0.25;
}
K[m * N + n] = S;
}
}
如果您关心S 的值,所有这些(当然)都很重要,因为有人可能会争辩说,由于您使用的函数会产生随机值,因此保持 S 值的顺序并不是最重要的.
带有线程安全随机生成器的版本
版本 1
#pragma omp parallel
{
unsigned int myseed = omp_get_thread_num();
#pragma omp for collapse(2)
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
for (int o = 0; o < O; o++) {
double r = ((double) rand_r(&myseed) / (double) RAND_MAX);
S += r - 0.25;
}
K[m * N + n] = S;
}
}
}
第 2 版
double *K = (double*) calloc(M * N, sizeof(double));
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
#pragma omp parallel
{
unsigned int myseed = omp_get_thread_num();
#pragma omp for reduction(+:S)
for (int o = 0; o < O; o++) {
double r = ((double) rand_r(&myseed) / (double) RAND_MAX);
S += r - 0.25;
}
}
K[m * N + n] = S;
}
}
编辑:
对原始问题进行更改。我想要平行和
顺序代码具有相同的结果。
而不是:
#pragma omp parallel for collapse(2)
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
#pragma omp for reduction(+:S)
for (int o = 0; o < O; o++) {
S += o;
}
K[m * N + n] = S;
}
}
做:
for (int m = 0; m < M; m++) {
for (int n = 0; n < N; n++) {
#pragma omp parallel for reduction(+:S)
for (int o = 0; o < O; o++) {
S += o;
}
K[m * N + n] = S;
}
}