c++ - 删除单个构造会导致执行不正确

Question

此代码按预期工作:

#include <iostream>
#include <cmath>
#include <omp.h>

//https://stackoverflow.com/questions/37970024/jacobi-relaxation-in-mpi
#define max(a, b) (a)>(b)?(a):(b)

const int m = 2001;
const int n = 1500;
const int p = 4;

double v[m + 2][m + 2];
double x[m + 2];
double y[m + 2];
double _new[m + 2][m + 2];
double maxdiffA[p + 1];
int icol, jrow;

int main() {
    omp_set_num_threads(p);

    double h = 1.0 / (n + 1);

    double start = omp_get_wtime();

    #pragma omp parallel for private(icol) shared(x, y, v, _new)
    for (icol = 0; icol <= n + 1; ++icol) {
        x[icol] = y[icol] = icol * h;

        _new[icol][0] = v[icol][0] = 6 - 2 * x[icol];

        _new[n + 1][icol] = v[n + 1][icol] = 4 - 2 * y[icol];

        _new[icol][n + 1] = v[icol][n + 1] = 3 - x[icol];

        _new[0][icol] = v[0][icol] = 6 - 3 * y[icol];
    }

    const double eps = 0.01;


    #pragma omp parallel private(icol, jrow) shared(_new, v, maxdiffA)
    {
        while (true) { //for [iters=1 to maxiters by 2]
            #pragma omp single
            for (int i = 0; i < p; i++) maxdiffA[i] = 0;
            #pragma omp for
            for (icol = 1; icol <= n; icol++)
                for (jrow = 1; jrow <= n; jrow++)
                    _new[icol][jrow] =
                            (v[icol - 1][jrow] + v[icol + 1][jrow] + v[icol][jrow - 1] + v[icol][jrow + 1]) / 4;
            #pragma omp for
            for (icol = 1; icol <= n; icol++)
                for (jrow = 1; jrow <= n; jrow++)
                    v[icol][jrow] = (_new[icol - 1][jrow] + _new[icol + 1][jrow] + _new[icol][jrow - 1] +
                                     _new[icol][jrow + 1]) / 4;

            #pragma omp for
            for (icol = 1; icol <= n; icol++)
                for (jrow = 1; jrow <= n; jrow++)
                    maxdiffA[omp_get_thread_num()] = max(maxdiffA[omp_get_thread_num()],
                                                         fabs(_new[icol][jrow] - v[icol][jrow]));

            #pragma omp barrier

            double maxdiff = 0.0;
            for (int k = 0; k < p; ++k) {
                maxdiff = max(maxdiff, maxdiffA[k]);
            }


            if (maxdiff < eps)
                break;
            #pragma omp single
            std::cout << maxdiff << std::endl;
        }
    }
    double end = omp_get_wtime();
    printf("start = %.16lf\nend = %.16lf\ndiff = %.16lf\n", start, end, end - start);

    return 0;
}

输出

1.12454
<106 iterations here>
0.0100436
start = 1527366091.3069999217987061
end = 1527366110.8169999122619629
diff = 19.5099999904632568

但是如果我删除

#pragma omp single
std::cout << maxdiff << std::endl;

程序要么似乎无限长地运行，要么我得到

start = 1527368219.8810000419616699
end = 1527368220.5710000991821289
diff = 0.6900000572204590

为什么会这样？

Answer 1

您在 while 循环开始时覆盖 maxdiffA - 这必须与在结束时读取 maxdiffA 以检查条件分开。否则，一个线程可能在另一个线程有机会读取它们之前就已经重置了这些值。由于 omp single 构造末尾的隐式屏障，循环末尾的 omp single 构造充当隔离。但是，omp single 结构的开头没有障碍。此外，“一大堆代码” 也不是安全屏障。因此，如果没有有效的隐式屏障，您必须使用 #pragma omp barrier 保护对重置代码的访问。

也就是说，我强烈建议重组代码以具有共享退出条件，该条件也在 single 构造中计算。这使得所有线程进程同时退出 while 循环变得更加清楚。否则代码定义错误。