dependencies - OpenMP : communicating values between threads 中的并行累积(前缀)总和

Question

假设我有一个函数 f(i)这取决于索引 i (以及其他无法预先计算的值)。
我想填充一个数组 a以便 a[n] = sum(f(i)) from i=0 to n-1 .

编辑:在 Hristo Iliev 发表评论后，我意识到我在做的是 cumulative/prefix sum .

这可以用代码编写为

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 )。

Answer 1

您可以将您的策略​​扩展到任意数量的子区域，并使用任务递归地减少它们:

#include<vector>
#include<iostream>

using namespace std;

const int n          = 10000;
const int baseLength = 100;

int f(int ii) {
  return ii;
}

int recursiveSumBody(int * begin, int * end){

  size_t length  = end - begin;
  size_t mid     = length/2;
  int    sum     = 0;


  if ( length < baseLength ) {
    for(size_t ii = 1; ii < length; ii++ ){
        begin[ii] += begin[ii-1];
    }
  } else {
#pragma omp task shared(sum)
    {
      sum = recursiveSumBody(begin    ,begin+mid);
    }
#pragma omp task
    {
      recursiveSumBody(begin+mid,end      );
    }
#pragma omp taskwait

#pragma omp parallel for
    for(size_t ii = mid; ii < length; ii++) {
      begin[ii] += sum;
    }

  }
  return begin[length-1];
}

void recursiveSum(int * begin, int * end){

#pragma omp single
  {
    recursiveSumBody(begin,end);
  }    
}


int main() {

  vector<int> a(n,0);

#pragma omp parallel
  {
    #pragma omp for
    for(int ii=0; ii < n; ii++) {          
      a[ii] = f(ii);
    }  

    recursiveSum(&a[0],&a[n]);

  }
  cout << n*(n-1)/2 << endl;
  cout << a[n-1] << endl;

  return 0;
}