Цикл OpenMP дает другой результат для одного и того же последовательного цикла

У меня есть серийный код:

double* a = malloc((1000000) * sizeof(double));
double* b = malloc((1000000) * sizeof(double));
double totalA = 0;

for (int i = 0; i < 1000000; i++) {

    if (i == 0) {
        a[i] = sin(i);
    }

    b[i] = sin(i+1);

    if (i < 1000000-1) {
        a[i+1] = b[i];
    }

    totalA += a[i];
}

Выход totalA после этого последовательного цикла равен 0.232883978073.

Затем у меня есть версия OpenMP (примечание: все переменные повторно инициализируются):

double* a = malloc((1000000) * sizeof(double));
double* b = malloc((1000000) * sizeof(double));
double totalA = 0;

#pragma omp parallel for reduction(+:totalA)
for (int i = 0; i < 1000000; i++) {

    if (i == 0) {
        a[i] = sin(i);
    }

    b[i] = sin(i+1);

    if (i < 1000000-1) {
        a[i+1] = b[i];
    }

    totalA += a[i];
}

Однако вывод totalA из этого кода равен -0.733714826779.

Я не могу понять для жизни меня, почему это отличается.

Спасибо.

ОБНОВЛЕНИЕ

После некоторой игры кажется, что операторы if в цикле странным образом игнорируются. Фактические операторы в блоках if выполняются на всех итерациях цикла (как если бы предложение if не существовало).

Например, изменение блока if на:

if (i < 555555) {
    a[i+1] = b[i];
}

кажется, что нет абсолютно никакой разницы.

Я до сих пор не знаю, что здесь происходит.


c openmp parallel-processing multicore reduction
person myles    schedule 28.11.2014    source источник
comment
Другой ответ набрал намного больше, чем мой. Пожалуйста, дайте ему принятый ответ.   -  person Z boson    schedule 03.12.2014

Ответы (2)


arrow_upward
5
arrow_downward

Ваш код содержит условие гонки. Конфликтующие операторы — это присваивание a[i+1] = b[i];, которое записывает в массив a, и оператор totalA += a[i];, который читает из a.

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

Чтобы еще больше продемонстрировать проблему, упорядочивание сегмента цикла, содержащего конфликтующие операторы, решает проблему (но, скорее всего, разрушительно для вашей производительности):

#pragma omp parallel for ordered reduction(+:totalA)    
for (int i = 0; i < 1000000; i++) {

   if (i == 0) {
     a[i] = sin(i);
   }

  b[i] = sin(i+1);

  #pragma omp ordered
  {
    if (i < 1000000-1) {
      a[i+1] = b[i];
    }

    totalA += a[i];
  }
}

Вероятно, лучше вообще избежать этой проблемы и переписать вашу программу, чтобы избавиться от зависимостей, переносимых циклом:

#define N 1000000

double* a = malloc(N * sizeof(double));
double* b = malloc(N * sizeof(double));
double totalA = 0;

a[0] = sin(0);
totalA += a[0];

#pragma omp parallel for reduction(+:totalA)
for (int i = 0; i < N - 1; i++) {
  b[i] = sin(i + 1);
  a[i + 1] = b[i];
  totalA += a[i + 1];
}

b[N - 1] = sin(N);

Наконец, обратите внимание, что, как и sin(0) ≈ 0.0, операторы

a[0] = sin(0);
totalA += a[0];

можно просто заменить на

a[0] = 0.0;
person Stefan Holdermans    schedule 28.11.2014
comment
Я не знаком с OpenMP и предполагал, что он сам об этом позаботится. - person 2501; 28.11.2014
comment
@ 2501, поместив #pragma omp for в цикл, вы утверждаете, что он может выполняться параллельно. Проверка безопасности этого утверждения является вашей работой, а не компилятора. (В общем случае компилятор не может этого сделать, так как в цикле могут быть вызовы отдельно скомпилированного кода, который, возможно, даже не был написан при запуске компилятора!) - person Jim Cownie; 01.12.2014
comment
Кажется, люди думают, что ваш ответ намного лучше моего (+5 — отличный ответ для OpenMP). Я ошибся в своем ответе? Мой ответ вводит в заблуждение? - person Z boson; 03.12.2014
comment
@Zboson Я думаю, твой ответ в порядке. Вы сразу же перешли к решению проблемы, где я больше сосредоточился на объяснении проблемы. Возможно, люди оценили это объяснение? Я бы сказал, что ответы хорошо дополняют друг друга. - person Stefan Holdermans; 03.12.2014
comment
Ответ Стефана объяснил мне, почему у меня возникла проблема (именно поэтому я проголосовал за нее). Однако ответ @Zboson «решил» мою проблему, поэтому я пометил его как «решение». Как вы упомянули Стефана, оба ответа дополняют друг друга и оба мне помогли. Спасибо вам обоим. - person myles; 03.12.2014
comment
Стефан, вы можете добавить простое исправление @Zboson в конец своего ответа, и я приму его как лучшее решение. Надеюсь, тогда у кого-нибудь еще, у кого есть эта проблема, будет как объяснение, так и решение? - person myles; 04.12.2014
comment
Я рад, что вы заменили totalA += b[i]; на totalA += a[i + 1];. Я хотел предложить это потом. Кроме того, я думаю, вы хотели написать sin(0) = 0.0 вместо ≈ 0.0. Или, может быть, вы имеете в виду sin(x) ≈ x вместо x << 1. - person Z boson; 04.12.2014

arrow_upward
3
arrow_downward

Вы можете упростить свой основной цикл до:

a[0] = sin(0);
for (int i = 0; i < N-1; i++) {
    b[i] = sin(i+1);
    a[i+1] = b[i];
    totalA += a[i];
}
totalA += a[N-1];
b[N-1] = sin(N);

Назовем итератор для thread1 i1 и thread2 i2. Затем, когда, например. i2=i1+1 запись в a[i1+1] и чтение в a[i2] будут иметь один и тот же адрес, а прочитанное или записанное значение будет зависеть от того, какой поток находится первым. Это состояние гонки.

Простое исправление состоит в том, чтобы заметить, что totalA = a[0] + b[0] + b[1] + ... b[N-2].

a[0] = sin(0);
totalA += a[0];
#pragma omp parallel for reduction(+:totalA)
for (int i = 0; i < N-1; i++) {
    b[i] = sin(i+1);
    a[i+1] = b[i];
    totalA += b[i];
}
b[N-1] = sin(N);

Это производит 0.232883978073 для N=1000000.

person Z boson    schedule 28.11.2014