OpenMP 并行减少会产生错误的结果

我正在使用信号矩阵，我的目标是计算一行所有元素的总和。该矩阵由以下结构表示：

typedef struct matrix {
  float *data;
  int rows;
  int cols;
  int leading_dim;
} matrix;

我不得不提到矩阵是按列优先顺序存储的（http://en.wikipedia.org/wiki/Row-major_order#Column-major_order http://en.wikipedia.org/wiki/Row-major_order#Column-major_order），这应该解释这个公式column * tan_hd.rows + row用于检索正确的索引。

for(int row = 0; row < tan_hd.rows; row++) {
    float sum = 0.0;
    #pragma omp parallel for reduction(+:sum)
    for(int column = 0; column < tan_hd.cols; column++) {
        sum += tan_hd.data[column * tan_hd.rows + row];
    }
    printf("row %d: %f", row, sum);
}

如果没有 OpenMP 编译指示，交付的结果是正确的，如下所示：

row 0: 8172539.500000 row 1: 8194582.000000 

一旦我添加#pragma omp...如上所述，返回不同（错误）的结果：

row 0: 8085544.000000 row 1: 8107186.000000

以我的理解，reduction(+:sum)创建私人副本sum对于每个线程，完成循环后，这些部分结果将被求和并写回全局变量sum再次。我做错了什么？

我很感激你的建议！

Use the 卡汉求和算法 http://en.wikipedia.org/wiki/Kahan_summation_algorithm

• 它具有与朴素求和相同的算法复杂性
• 它将大大提高求和的准确性，而无需将数据类型切换为双倍。

通过重写代码来实现它：

for(int row = 0; row < tan_hd.rows; row++) {
    float sum = 0.0, c = 0.0;
    #pragma omp parallel for reduction(+:sum, +:c)
    for(int column = 0; column < tan_hd.cols; column++) {
        float y = tan_hd.data[column * tan_hd.rows + row] - c;
        float t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    sum = sum - c;
    printf("row %d: %f", row, sum);
}

您还可以切换所有float to double达到更高的精度，但由于你的数组是float数组中，应该只存在末尾有效数字的数量差异。