c - 用于减少的 OpenMP 并行提供了错误的结果

Question

我正在处理一个信号矩阵，我的目标是计算一行中所有元素的总和。矩阵由以下结构表示:

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

我必须提到该矩阵以列优先顺序 ( 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 pragma，交付的结果是正确的，如下所示:

row 0: 8172539.500000 row 1: 8194582.000000 

一旦我如上所述添加 #pragma omp...，就会返回一个不同的(错误的)结果:

row 0: 8085544.000000 row 1: 8107186.000000

在我的理解中，reduction(+:sum) 为每个线程创建了 sum 的私有(private)副本，并在完成循环后将这些部分结果相加并写回再次到全局变量 sum。我做错了什么？

感谢您的建议!

Answer 1

使用 Kahan summation algorithm

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

通过重写代码来实现它:

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 切换为 double 以获得更高的精度，但是由于您的数组是 float 数组，因此应该只有是末尾有效数字数量的差异。