Этот фрагмент кода решает систему линейных уравнений методом простой итерации:
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) сложно реализовать динамическую балансировку
Простой
