Как умно распараллелить вложенный цикл OpenMP?

Question

[aab137]

Этот фрагмент кода решает систему линейных уравнений методом простой итерации:

int n_iter;
    double error = eps + 1;
    double *x_new = (double *)malloc(N * sizeof(double));
    double *Ax = (double *)malloc(N * sizeof(double));
    for (n_iter = 0; error >= eps; n_iter++) {
        error = 0;
        init_vector(Ax, N, 0);

        for (int i = 0; i < N; i++) {  // находим вектор Ax
            for (int j = 0; j < N; j++) {
                Ax[i] += A[i * N + j] * x[j];
            }
        }

        for (int i = 0; i < N; i++) {  // находим новый x и значение ошибки
            x_new[i] = x[i] - tau * (Ax[i] - b[i]);
            error += pow(x_new[i] - x[i], 2);
        }
        error = sqrt(error);
        
        // Обновление вектора x
        double *t = x;
        x = x_new;
        x_new = t;
    }

Нужно его распараллелить, создав одну параллельную секцию `#pragma omp parallel`, охватывающую весь итерационный алгоритм.
Почему просто не распараллелить два вложенных цикла - потому что тогда разбиение программы на процессы будет отнимать время на каждой итерации внешнего цикла.
Как это сделать - не понимаю. Нужно, чтобы внешний цикл сохранил порядок выполнения итераций, а вложенные циклы были распараллелены, при этом должна быть одна `#pragma omp parallel` перед внешним циклом.
Есть возможность разбить итерации по процессам вручную, как в этом примере:

#pragma omp parallel private(thread_id)
    {
        thread_id = omp_get_thread_num();
        for (iter_count = 0; accuracy > EPSILON && iter_count < MAX_ITERATION_COUNT; ++iter_count)
        {
            calc_Axb(A + line_offsets[thread_id] * N, x, b + line_offsets[thread_id], 
                     Axb + line_offsets[thread_id], line_counts[thread_id]);
#pragma omp barrier

            calc_next_x(Axb + line_offsets[thread_id], x + line_offsets[thread_id], 
                        TAU, line_counts[thread_id]);

#pragma omp single
            accuracy = 0;

#pragma omp atomic
            accuracy += calc_norm_square(Axb + line_offsets[thread_id], line_counts[thread_id]);
#pragma omp barrier

#pragma omp single
            accuracy = sqrt(accuracy) / b_norm;
        }
    }

Но этот способ мне не нравится, т.к. 1) выглядит костыльным, 2) сложно реализовать динамическую балансировку

Answer 1

[Wataru]

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

Внутри можно циклы по i объединить все в один.

Ну и параллельте цикл по i через pragma omp for.

Comments

[aab137]

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

[Wataru]

aab137, Ну хорошо, помеяйте местами циклы по i и по n_iter. Но надо будет добавить добавить barrier в конце каждой итерации цикла по n_iter. Так создадуться куча потоков, но они будут ждать друг друга каждую итерацию.