Удаление одной конструкции приводит к некорректному выполнению

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

Вы перезаписываете maxdiffA в начале цикла while - это должно быть изолировано от чтения maxdiffA в конце, чтобы проверить условие.В противном случае один поток может уже сбросить значения, прежде чем другой поток получит возможность прочитать их.Конструкция omp single в конце цикла действует как изоляция из-за неявного барьера в конце конструкции omp single.Однако в начале конструкции omp single барьер отсутствует.Также «много кода» не является безопасным барьером.Поэтому, если нет действительного неявного барьера, вы должны защитить запись в коде сброса с помощью #pragma omp barrier.

При этом я настоятельно рекомендую реструктурировать код, чтобы иметь общее условие выхода, которое также вычисляется вsingle конструкция.Это делает более ясным, что все процессы потока выходят из цикла while одновременно.В противном случае код не определен.