Как распараллелить эволюцию 2D времени

У меня есть код C ++, который выполняет эволюцию во времени четырех переменных, которые живут в двумерной пространственной сетке. Чтобы сэкономить время, я попытался распараллелить мой код с OpenMP, но я просто не могу заставить его работать: независимо от того, сколько ядер я использую, время выполнения остается в основном одинаковым или увеличивается. (Мой код использует 24 ядра или столько, сколько я укажу, так что компиляция не проблема.)

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

Макет моего кода:

for (int t = 0; t < max_time_steps; t++) {

// do some book-keeping
...

// perform time step
// (1) calculate righthand-side of ODE:
for (int i = 0; i < nr; i++) {
for (int j = 0; j < ntheta; j++) {
rhs[0][i][j] = A0[i][j] + B0[i][j] + ...;
rhs[1][i][j] = A1[i][j] + B1[i][j] + ...;
rhs[2][i][j] = A2[i][j] + B2[i][j] + ...;
rhs[3][i][j] = A3[i][j] + B3[i][j] + ...;
}
}

// (2) perform Euler step (or Runge-Kutta, ...)
for (int d = 0; d < 4; d++) {
for (int i = 0; i < nr; i++) {
for (int j = 0; j < ntheta; j++) {
next[d][i][j] = current[d][i][j] + time_step * rhs[d][i][j];
}
}
}

}

Я думал, что этот код должен быть довольно легко распараллелить … Я поставил «#pragma omp parellel for» перед циклами (1) и (2), и я также указал количество ядер (например, 4 ядра для цикла ( 2) поскольку существует четыре переменные), но ускорения просто нет.

Я обнаружил, что OpenMP довольно умно относится к тому, когда создавать / уничтожать потоки. То есть он понимает, что скоро снова понадобятся потоки, и тогда они только засыпают, чтобы сэкономить время.

Я думаю, что одна «проблема» заключается в том, что мой временной шаг закодирован в подпрограмме (я использую RK4 вместо Эйлера), а вычисление правой части снова в другой подпрограмме, которая вызывается функцией time_step (). Итак, я считаю, что из-за этого OpenMP не может видеть, что потоки должны оставаться открытыми дольше, и, следовательно, потоки создаются и уничтожаются на каждом временном шаге.

Было бы полезно поместить «#pragma omp parallel» перед циклом времени, чтобы потоки создавались в самом начале? А затем происходит ли фактическое распараллеливание для правой части (1) и шага Эйлера (2)? Но как мне это сделать?

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

Я сейчас удалил d циклы (делая индексы явными) и свернул i а также j циклы (при запуске по всему 2D-массиву только с одной переменной).

Код выглядит так:

for (int t = 0; t < max_time_steps; t++) {

// do some book-keeping
...

// perform time step
// (1) calculate righthand-side of ODE:
#pragma omp parallel for
for (int i = 0; i < nr*ntheta; i++) {
rhs[0][0][i] = A0[0][i] + B0[0][i] + ...;
rhs[1][0][i] = A1[0][i] + B1[0][i] + ...;
rhs[2][0][i] = A2[0][i] + B2[0][i] + ...;
rhs[3][0][i] = A3[0][i] + B3[0][i] + ...;
}

// (2) perform Euler step (or Runge-Kutta, ...)
#pragma omp parallel for
for (int i = 0; i < nr*ntheta; i++) {
next[0][0][i] = current[0][0][i] + time_step * rhs[0][0][i];
next[1][0][i] = current[1][0][i] + time_step * rhs[1][0][i];
next[2][0][i] = current[2][0][i] + time_step * rhs[2][0][i];
next[3][0][i] = current[3][0][i] + time_step * rhs[3][0][i];
}

}

Размер nr*ntheta является 400*40=1600 и я сделаю max_time_steps=1000 временные шаги. Тем не менее, распараллеливание не приводит к ускорению:

Время выполнения без OpenMP (результат time в командной строке):

real   0m23.597s
user   0m23.496s
sys    0m0.076s

Время выполнения с OpenMP (24 ядра)

real   0m23.162s
user   7m47.026s
sys    0m0.905s

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

Одна особенность, которую я не показываю в приведенном выше фрагменте кода, состоит в том, что мои переменные на самом деле не являются doubleс, но самоопределенная структура из двух doubleс, которые напоминают реальную и мнимую части. Но я думаю, что это не должно иметь значения.