Supongamos que tengo una funciónf(i) lo que depende de un índicei (entre otros valores que no pueden ser precomputados). Quiero llenar una matriza así que esoa[n] = sum(f(i)) from i=0 to n-1.

Editar: Después de un comentario de Hristo Iliev, me di cuenta de que lo que estoy haciendo es unsuma acumulativa / prefijo.

Esto se puede escribir en código como

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

Ahora quiero usar OpenMP para hacer esto en paralelo. Una forma en que podría hacer esto con OpenMP es escribir los valores paraf(i) En paralelo y luego cuidar la dependencia en serie. Sif(i) es una función lenta, entonces esto podría funcionar bien ya que el bucle no paralelo es simple.

#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];
}

Pero es posible hacer esto sin el bucle no paralelo con OpenMP. La solución, sin embargo, que he encontrado es complicada y tal vez un poco complicada. Entonces, mi pregunta es si hay una forma más sencilla y complicada de hacer esto con OpenMP.

El siguiente código básicamente ejecuta el primer código que enumeré para cada hilo. El resultado es que los valores dea En un hilo dado son correctos hasta una constante. Guardo la suma de cada hilo en una matrizsuma connthreads+1 elementos. Esto me permite comunicarme entre hilos y determinar el desplazamiento constante para cada hilo. Entonces corrijo los valores dea[i] con el offset.

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;

Una prueba simple es sólo para establecerf(i) = i. Entonces la solución esa[i] = i*(i+1)/2 (y en el infinito es-1/12).