OpenMP并行for循环归约得出错误结果。

Question

OpenMP并行for循环归约得出错误结果。

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。那我做错了什么呢？

感谢您的建议！

- Max Plauth

我认为原因是column是局部变量。对于每个并行操作，column都会被初始化为零，但实际上不应该这样做。将其移出for循环即可。 - MYMNeo

浮点数加法不是结合律——如果更改元素求和的顺序，则会以不同的方式累积小的舍入误差。在串行情况下可能使用64位或80位内部精度，但在规约阶段会丢失每个值被转换为单精度的精度。 - Hristo Iliev

@Hristo Iliev：谢谢，我不知道浮点数不是可结合的。只是为了好玩，我会尝试使用双精度浮点数，观察误差是否变小或保持不变 :) - Max Plauth

第一次尝试使用double值只会导致更大的误差。另一方面，在omp部分期间，tan_hd.data和row都没有被更改/写入，omp应该负责column... - Max Plauth

@HighPerformanceMark，我也没有看到任何数据竞争的迹象，这就是为什么我猜测可能是由于舍入误差引起的。拥有许多正负值的数组很容易导致灾难性的取消效应。 - Hristo Iliev

显示剩余2条评论

1个回答

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Kyle_the_hacker · Accepted Answer

使用Kahan求和算法

它具有与朴素求和相同的算法复杂度
它将极大地增加求和的精度，而无需切换数据类型为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 数组，在结尾处只会有显著数字的数量上的差异。