Почему OpenMP не может суммировать эти цифры?

Question

Рассмотрим следующий минимальный пример кода C. При компиляции и выполнении с export OMP_NUM_THREADS=4 && gcc -fopenmp minimal2.c && ./a.out (homebrew GCC 5.2.0 на OS X 10.11) это обычно приводит к правильному поведению, т. Е. Семи строк с одинаковым числом. Но иногда это случается:

[ ] 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 <stdio.h> 
#include <stdlib.h> 
#include <time.h> 

#ifdef _OPENMP 
#include <omp.h> 
#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; 
} 

EDIT: Жиль обратил мое внимание, что ошибки, начиная с 15-й значащей цифры являются нормальными, как я переоценил точность двойной. Я также не могу воспроизвести ошибочное поведение с 2x правильным номером на машине Debian, так что это может быть родной gcc или Mac.

У меня была проблема с аналогичным выпуском here, но эти два, кажется, не связаны (по крайней мере, в моих глазах), поэтому я начал это как отдельный вопрос.

Answer 1

Я сильно подозреваю, что это потому, что floating-point addition is not associative. В результате OpenMP суммирует умножения в разных порядках, давая несколько разные результаты.

OpenMP 4.0 spec, section 1.3 Execution Model говорит:

Например, последовательное снижение добавления может иметь различную структуру аддитивных ассоциаций, чем параллельное сокращение. Эти различные ассоциации могут изменять результаты сложения с плавающей запятой.

См. OpenMP parallel for reduction delivers wrong results за предлагаемое решение.