c - OpenMP 并行 for 循环问题与指向指针的共享指针

Question

我刚刚学习 OpenMP，在并行化一个 for 循环时遇到了问题，我很确定它符合并行化的标准。下面的 SSCCE 展示了我在实际代码中遇到的问题。

问题是，当我用 OpenMP 编译和运行它时，wr[i][j] 和 wi[i][j] 的内容总是 1.00和 0.00 分别。如果我注释掉 pragma 并连续运行代码，则 wr 和 wi 包含我期望的值。即使我在用 OpenMP 编译后只用 1 个线程运行它，代码仍然会出现问题。必须注释掉 pragma 才能连续运行它。

我认为我的问题与 wr 和 wi 是指向指针的指针有关。如果是这种情况，我该如何解决？如果那不是问题是什么？非常感谢任何见解。

// compilation gcc -fopenmp -o test test.c -lm 
// usage ./test <Problem Size (Must be a power of 2)> <Num Threads>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

int main (int argc, char *argv[]){
    double **wr, **wi;  
    double *data;
    double wtemp, wpr, wpi, theta;
    int i, j, n, nn, mmax, istep, num_mmax_iterations, T;
    const int isign = -1;

    if(argc < 3){
        printf("ERROR, usage: %s <N (must be a power of 2)> <T>\n", argv[0]);
        exit(1);
    } else {
        n = atoi(argv[1]);
        T = atoi(argv[2]);
    }

    num_mmax_iterations = log2(n);
    wr = malloc(sizeof(double *) *num_mmax_iterations);
    wi = malloc(sizeof(double *) *num_mmax_iterations);

    mmax = 2;
    for(i = 0; i < num_mmax_iterations; i++) {
        wr[i] = malloc(sizeof(double) * (mmax >> 1));
        wi[i] = malloc(sizeof(double) * (mmax >> 1));
        wr[i][0] = 1.0;
        wi[i][0] = 0.0;
        mmax <<= 1;
    }

    mmax = 2;
    #pragma omp parallel for private(i, j, theta, wtemp, wpr, wpi) firstprivate(mmax) shared(num_mmax_iterations, wr, wi) num_threads(T)
    for(i = 0; i < num_mmax_iterations; i++) {
        theta = isign * (6.28318530717959 / mmax);
        wtemp = sin(0.5 * theta);
        wpr = -2.0 * wtemp * wtemp;
        wpi = sin(theta);       

        for(j = 1; j < (mmax >> 1); j++) {
            wr[i][j] = wr[i][j - 1] * wpr - wi[i][j - 1] * wpi + wr[i][j - 1];
            wi[i][j] = wi[i][j - 1] * wpr + wr[i][j - 1] * wpi + wi[i][j - 1];
        }

        mmax <<= 1;
    }

    for(i = 0; i < num_mmax_iterations; i++) {      
        for(j = 0; j < (mmax >> 1); j++) {
            printf("%.2f | %.2f\n", wr[i][j], wi[i][j]);
        }
    }

    for(i = 0; i < num_mmax_iterations; i++) {
        free(wr[i]);
        free(wi[i]);
    }
    free(wr);
    free(wi);
    return 0;
}

编辑

运行 ./test 8 1 时删除编译指示的预期输出

1.00 | 0.00
0.00 | 0.00
0.00 | 1.00
0.00 | -0.00
1.00 | 0.00
-0.00 | -1.00
0.00 | 1.00
-1.00 | 0.71
1.00 | -0.00
0.71 | -0.71
-0.00 | 0.00
-0.71 | -0.71
0.00 | -1.00
-0.71 | -0.71
-1.00 | 1.00
-0.71 | 0.92
1.00 | 0.00
-0.00 | -1.00
0.00 | 1.00
-1.00 | 0.71
1.00 | -0.00
0.71 | -0.71
-0.00 | 0.00
-0.71 | -0.71
0.00 | -1.00
-0.71 | -0.71
-1.00 | 1.00
-0.71 | 0.92
1.00 | 0.71
0.92 | 0.38
0.71 | -0.00
0.38 | -0.38
1.00 | 0.00
0.71 | -0.71
-0.00 | -1.00
-0.71 | -0.71
0.00 | 1.00
-0.71 | 0.92
-1.00 | 0.71
-0.71 | 0.38
1.00 | -0.00
0.92 | -0.38
0.71 | -0.71
0.38 | -0.92
-0.00 | 0.00
-0.38 | -0.38
-0.71 | -0.71
-0.92 | -0.92
1.00 | 0.00
0.92 | -0.38
0.71 | -0.71
0.38 | -0.92
-0.00 | -1.00
-0.38 | -0.92
-0.71 | -0.71
-0.92 | -0.38
0.00 | 0.00
-0.38 | 0.00
-0.71 | 0.00
-0.92 | 0.00
-1.00 | 0.00
-0.92 | 0.00
-0.71 | 0.00
-0.38 | 0.00

运行 ./test 8(任意数量的线程)时完整保留 pragma 的实际输出

1.00 | 0.00
1.00 | 0.00
1.00 | 0.00

注意:实际问题的规模要大得多，预期数据的变化通常比上面的大。

Answer 1

当代码连续运行时，值 mmax 在外循环的每次迭代中都会被修改。您的代码需要修改，以便每个线程为自己计算正确的 mmax 值。此外，外部循环变量 i 应该在线程之间共享，而不是私有(private)的。这是对您的 openMP 循环的建议，无论是否使用 pragma，它都会给我相同的结果。

#pragma omp parallel for private( j, theta, wtemp, wpr, wpi,mmax) num_threads(T)
for(i = 0; i < num_mmax_iterations; i++) {
    mmax = pow(2,i+1);
    theta = isign * (6.28318530717959 / mmax);
    wtemp = sin(0.5 * theta);
    wpr = -2.0 * wtemp * wtemp;
    wpi = sin(theta);       

    for(j = 1; j < (mmax >> 1); j++) {
        wr[i][j] = wr[i][j - 1] * wpr - wi[i][j - 1] * wpi + wr[i][j - 1];
        wi[i][j] = wi[i][j - 1] * wpr + wr[i][j - 1] * wpi + wi[i][j - 1];
    }

}
// mmax has to be set at the end since there is no guarantee of which thread will finish last
mmax = pow(2,1+num_mmax_iterations);