title: Parallel cumulative (prefix) sums in OpenMP: communicating values between threads
body: <p>Assume I have a function <code>f(i)</code> which depends on an index <code>i</code> (among other values which cannot be precomputed). 
  I want to fill an array <code>a</code> so that <code>a[n] = sum(f(i)) from i=0 to n-1</code>.</p>
  
  <p><strong>Edit:</strong> After a comment by Hristo Iliev I realized what I am doing is a <a href="http://en.wikipedia.org/wiki/Prefix_sum" rel="nofollow">cumulative/prefix sum</a>.</p>
  
  <p>This can be written in code as</p>
  
  <pre><code>float sum = 0;
  for(int i=0; i&lt;N; i++) {
      sum += f(i);
      a[i] = sum;
  }
  </code></pre>
  
  <p>Now I want to use OpenMP to do this in parallel.  One way I could do this with OpenMP is to write out the values for <code>f(i)</code> in parallel and then take care of the dependency in serial.  If <code>f(i)</code> is a slow function then this could work well since the non-paralleled loop is simple.</p>
  
  <pre><code>#pragma omp parallel for
  for(int i=0; i&lt;N; i++) {
      a[i] = f(i);
  }
  for(int i=1; i&lt;N; i++) {
      a[i] += a[i-1];
  }
  </code></pre>
  
  <p>But it's possible to do this without the non-parallel loop with OpenMP.  The solution, however, that I have come up with is complicated and perhaps hackish.  So my question is if there is a simpler less convoluted way to do this with OpenMP?  </p>
  
  <p>The code below basically runs the first code I listed for each thread.  The result is that values of <code>a</code> in a given thread are correct up to a constant.  I save the sum for each thread to an array <code>suma</code> with <code>nthreads+1</code> elements.  This allows me to communicate between threads and determine the constant offset for each thread.  Then I correct the values of <code>a[i]</code> with the offset.</p>
  
  <pre><code>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&lt;finish; i++) {
          sum += f(i);
          a[i] = sum;
      }
      suma[ithread+1] = sum;
      #pragma omp barrier
      float offset = 0;
      for(int i=0; i&lt;(ithread+1); i++) {
          offset += suma[i];
      }
      for(int i=start; i&lt;finish; i++) {
          a[i] += offset;
      }
  }
  delete[] suma;
  </code></pre>
  
  <p>A simple test is just to set <code>f(i) = i</code>.  Then the solution is <code>a[i] = i*(i+1)/2</code> (and at infinity it's <a href="http://en.wikipedia.org/wiki/Ramanujan_summation#Sum_of_divergent_series" rel="nofollow">-1/12</a>).</p>