Параллельные кумулятивные (префиксные) суммы в OpenMP: передача значений между потоками

Предположим, у меня есть функцияf(i) который зависит от индексаi (среди других значений, которые не могут быть предварительно вычислены). Я хочу заполнить массивa так чтоa[n] = sum(f(i)) from i=0 to n-1.

Редактировать: После комментария Христо Илиева я понял, что я делаю, этонакопленная / префиксная сумма.

Это можно записать в коде как

float sum = 0;
for(int i=0; i<N; i++) {
    sum += f(i);
    a[i] = sum;
}

Теперь я хочу использовать OpenMP, чтобы сделать это параллельно. Один из способов сделать это с помощью OpenMP - выписать значения дляf(i) параллельно, а затем позаботиться о зависимости в последовательном. Еслиf(i) это медленная функция, тогда это может хорошо работать, так как непараллельный цикл прост.

#pragma omp parallel for
for(int i=0; i<N; i++) {
    a[i] = f(i);
}
for(int i=1; i<N; i++) {
    a[i] += a[i-1];
}

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

Код ниже в основном запускает первый код, который я перечислил для каждого потока. Результатом является то, что значенияa в заданном потоке верны с точностью до константы. Я сохраняю сумму для каждого потока в массивsuma с участиемnthreads+1 элементы. Это позволяет мне общаться между потоками и определять постоянное смещение для каждого потока. Тогда я исправляю значенияa[i] со смещением.

float *suma;
#pragma omp parallel
{
    const int ithread = omp_get_thread_num();
    const int nthreads = omp_get_num_threads();
    const int start = ithread*N/nthreads;
    const int finish = (ithread+1)*N/nthreads;
    #pragma omp single
    {
        suma = new float[nthreads+1];
        suma[0] = 0;
    }
    float sum = 0;
    for (int i=start; i<finish; i++) {
        sum += f(i);
        a[i] = sum;
    }
    suma[ithread+1] = sum;
    #pragma omp barrier
    float offset = 0;
    for(int i=0; i<(ithread+1); i++) {
        offset += suma[i];
    }
    for(int i=start; i<finish; i++) {
        a[i] += offset;
    }
}
delete[] suma;

Простой тест просто установитьf(i) = i, Тогда решениеa[i] = i*(i+1)/2 (и в бесконечности это-1/12).

Ответы на вопрос(1)

Ваш ответ на вопрос