OpenMP 程序执行中的错误答案

【问题标题】：error in OpenMP program executionOpenMP 程序执行中的错误
【发布时间】：2016-06-11 01:33:50
【问题描述】：

我正在尝试对在同一代码中初始化的数组的所有元素求和。由于每个元素都是相互独立的，因此我尝试并行执行求和。我的代码如下所示：

int main(int argc, char** argv) 
{
  cout.precision(20);
  double sumre=0.,Mre[11];
  int n=11;
  for(int i=0; i<n; i++)
   Mre[i]=2.*exp(-10*M_PI*i/(1.*n));

  #pragma omp parallel for reduction(+:sumre)
  for(int i=0; i<n; i++)
  {
   sumre+=Mre[i];
  }
  cout<<sumre<<"\n";
}

我编译和运行的：

g++ -O3 -o sum sumparallel.cpp -fopenmp
./sum

分别。我的问题是每次运行时输出都不同。有时它给了

2.1220129388411006488

或

2.1220129388411002047 有谁知道这里发生了什么？

【问题讨论】：

double 只有 15–17 位有效十进制数字精度。
我主要担心的是我对获得良好的精度非常感兴趣，所以我不知道如何使用 OpenMP 来管理它。
在几个 SO Q 和 As 中也相关并涵盖了 f-p 算术的非结合性和并行 f-p 算术的非确定性。我自己对知识总和的贡献之一——stackoverflow.com/questions/32074644/…——还有其他人。

标签： openmp

【解决方案1】：

其中一些 cmets 暗示了这里的问题，但可能存在两个不同的问题

双精度数没有 20 位小数精度

如果要打印sumre 的最大精度，请使用类似这样的内容

#include <float.h>
int maint(int argc, char* argv[])
{
  ...

  printf("%.*g", DBL_DECIMAL_DIG, number);

  return 0;
}

浮点运算是不可交换的

此属性的影响是舍入误差。事实上，您定义的函数（高斯函数）特别容易舍入求和。考虑到 OpenMP parallel for 的工作负载分布是未定义的，您在运行它时可能会收到不同的答案。要解决此问题，您可以使用kahan summation algorithm。使用 OpenMP 的实现如下所示：

  ...
  double sum = 0.0, c = 0.0; 

  #pragma omp parallel for reduction(+:sum, +:c) 
  for(i = 0; i < n; i++) 
  { 
    double y = Mre[i] - c; 
    double t = sum + y; 
    c = (t - sum) - y; 
    sum = t; 
  } 

  sum = sum - c;  
  ...

【讨论】：