为什么OpenMP无法对这些数字求和?

考虑以下最小C代码示例。 在使用export OMP_NUM_THREADS=4 && gcc -fopenmp minimal2.c && ./a.out X 10.11上的自制GCC 5.2.0)进行编译和执行时,这通常会产生正确的行为,即七行具有相同的编号。 但有时会发生这种情况:

 [ ] bsum=1.893293142303100e+03 [1] asum=1.893293142303100e+03 [2] asum=1.893293142303100e+03 [0] asum=1.893293142303100e+03 [3] asum=3.786586284606200e+03 [ ] bsum=1.893293142303100e+03 [ ] asum=3.786586284606200e+03 equal: 0 

它看起来像一个竞争条件,但我的代码对我来说似乎很好。 我究竟做错了什么?

 #include  #include  #include  #ifdef _OPENMP #include  #define ID omp_get_thread_num() #else #define ID 0 #endif #define N 1400 double a[N]; double verify() { int i; double bsum = 0.0; for (i = 0; i < N; i++) { bsum += a[i] * a[i]; } fprintf(stderr, "[ ] bsum=%.15e\n", bsum); return bsum; } int main(int argc, char *argv[]) { int i; double asum = 0.0, bsum; srand((unsigned int)time(NULL)); //srand(1445167001); // fails on my machine for (i = 0; i < N; i++) { a[i] = 2 * (double)rand()/(double)RAND_MAX; } bsum = verify(); #pragma omp parallel shared(asum) { #pragma omp for reduction(+: asum) for (i = 0; i < N; i++) { asum += a[i] * a[i]; } fprintf(stderr, "[%d] asum=%.15e\n", ID, asum); } bsum = verify(); fprintf(stderr, "[ ] asum=%.15e\n", asum); return 0; } 

编辑:吉尔斯引起我的注意,从第15位有效数字开始的错误是正常的,因为我高估了双精度。 我也无法在Debian机器上以2x正确的数字重现错误行为,因此这可能是自制的gcc或Mac相关。

我在这里遇到类似问题的问题,但两者似乎没有关系(至少在我看来),所以我把它作为一个单独的问题开始。

我强烈怀疑这是因为浮点加法不是关联的 。 因此,OpenMP将不同顺序的乘法相加,得到的结果略有不同。

OpenMP 4.0规范1.3节执行模型说:

例如,串行加法减少可以具有与并行减少不同的加法关联模式。 这些不同的关联可能会改变浮点加法的结果。

请参阅OpenMP parallel for reduction为建议的解决方案提供错误的结果 。