3

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

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 编译指示，传递的结果是正确的，如下所示：

row 0: 8172539.500000 row 1: 8194582.000000

一旦我添加了#pragma omp...上述内容，就会返回一个不同的（错误的）结果：

row 0: 8085544.000000 row 1: 8107186.000000

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

我很感激你的建议！

4

1 回答 1

2

使用Kahan 求和算法

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

通过重写代码来实现它：

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数组，所以最后应该只有有效数字的数量有所不同。

于 2013-08-02T12:24:44.070 回答