OpenMP并行化GEMM的错误结果

Question

我知道OpenMP在所有工作人员之间共享在外部作用域中声明的所有变量。这就是我的问题的答案。但我真的不明白为什么函数omp3提供正确的结果，而函数omp2提供错误的结果。

void omp2(double *A, double *B, double *C, int m, int k, int n) {
    for (int i = 0; i < m; ++i) {
#pragma omp parallel for
        for (int ki = 0; ki < k; ++ki) {
            for (int j = 0; j < n; ++j) {
                C[i * n + j] += A[i * k + ki] * B[ki * n + j];
            }
        }
    }
}

void omp3(double *A, double *B, double *C, int m, int k, int n) {
    for (int i = 0; i < m; ++i) {
        for (int ki = 0; ki < k; ++ki) {
#pragma omp parallel for
            for (int j = 0; j < n; ++j) {
                C[i * n + j] += A[i * k + ki] * B[ki * n + j];
            }
        }
    }
}

Answer 1

问题是，这条线有一个比赛条件：

C[i * n + j] += ...

不同的线程可以同时读写相同的内存位置(C[i * n + j])，从而导致数据竞争。在omp2中，这种数据竞争可能发生，但在omp3中不会发生。

解决方案是(正如@Victor Eijkhout所建议的)：重新排序您的循环，使用局部变量计算最内部循环的和。在本例中，C[i * n + j]只更新一次，因此您摆脱了数据竞争，最外层的循环可以并行化(这提供了最好的性能)：

   #pragma omp parallel for
    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            double sum=0;            
            for (int ki = 0; ki < k; ++ki) {            
                sum += A[i * k + ki] * B[ki * n + j];
            }
            C[i * n + j] +=sum;
        }        
    }

请注意，您可以使用collapse(2)子句，这可能会提高性能。