Как openmp распараллелить для цикла, который увеличивает две переменные

Question

В качестве первого шага в OpenMP я поставил перед собой задачу распараллелить некоторый алгоритм декомпозиции матрицы.Я выбрал Crout с поворотом, источник можно найти здесь: http://www.mymathlib.com/c_source/matrices/linearsystems/crout_pivot.c

В нижней части этой функции разложения есть внешний цикл for, который проходит по i и p_row одновременно.Конечно, OpenMP так же запутан, как и я, когда смотрит на это, и отказывается что-либо с ним делать.

После того, как я обдумал его, я думаю, что я распутал его в читабельную форму:

p_row = p_k + n;
for (i = k+1; i < n; i++) {
    for (j = k+1; j < n; j++) *(p_row + j) -= *(p_row + k) * *(p_k + j);
    p_row += n;
}

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

p_row = p_k + n;
#pragma omp parallel for private (i,j) shared (n,k,p_row,p_k)
for (i = k+1; i < n; i++) {
    for (j = k+1; j < n; j++) *(p_row + j) -= *(p_row + k) * *(p_k + j);
#pragma omp critical
    p_row += n;
#pragma omp flush(p_row)
}

И все же результаты в основном случайные.

Чего мне не хватает?

Answer 1

Я не проверял вашу адаптацию исходного кода, но у вашей программы есть несколько проблем.

#pragma omp parallel for private (i,j) shared (n,k,p_row,p_k)
Поведение по умолчанию - объявить переменные вне области общего доступа, поэтому объявление общего доступабесполезный.Но эти переменные должны , а не быть общими и отображаться как приватные.
n не изменяется во время итераций, поэтому лучше иметь локальную копию
то же для k и p_k
p_row, но вы очень хочется несколько копий p_row.Это то, что обеспечит правильную параллельную обработку, так что каждый поток обрабатывает разные строки.Проблема состоит в том, чтобы вычислить значение p_row в разных потоках.

Во внешнем цикле итерация 0 будет использовать p_row, вторая итерация p_row + n, итерация l p_row + l * n.Ваши итерации будут распределены по нескольким потокам.Предположим, что каждый поток обрабатывает m итераций.Поток 0 будет обрабатывать i = k + 1 до i = m + (k + 1) и от p_row до p_row + m * n, поток 1 i = m + 1 + (k + 1) до i = 2m + (k + 1) иот p_row + n * (m + 1) до p_row + 2 * m * n и т. д. Следовательно, вы можете вычислить значение, которое должно иметь p_row в начале цикла, со значением i.Вот возможная реализация

p_row = p_k + n;
#pragma omp parallel for private(i,j) firstprivate(n, k, p_row, p_k)
// first private insures initial values are kept
{ 
  for (i = k+1, p_row=p_row+(i-(k+1))*n; i < n; i++, p_row += n) {
    for (j = k+1; j < n; j++) 
       *(p_row + j) -= *(p_row + k) * *(p_k + j);
}

Инкремент p_row находится в цикле for.Это должно продолжать работать в последовательной среде.

Критика бесполезна (и в предыдущем коде была ошибка).Flush подразумевается в конце параллельного раздела (а прагма просто "omp flush").