title: OpenMP/C++: Parallel for loop with reduction afterwards - best practice?
body: <p>Given the following code...</p>
  
  <pre><code>for (size_t i = 0; i &lt; clusters.size(); ++i)
  {
      const std::set&lt;int&gt;&amp; cluster = clusters[i];
      // ... expensive calculations ...
      for (int j : cluster)
          velocity[j] += f(j); 
  } 
  </code></pre>
  
  <p>...which I would like to run on multiple CPUs/cores. The function <code>f</code> does not use <code>velocity</code>.</p>
  
  <p>A simple <code>#pragma omp parallel for</code> before the first for loop will produce unpredictable/wrong results, because the <code>std::vector&lt;T&gt; velocity</code> is modified in the inner loop. Multiple threads may access and (try to) modify the same element of <code>velocity</code> at the same time.</p>
  
  <p>I think the first solution would be to write <code>#pragma omp atomic</code> before the <code>velocity[j] += f(j);</code>operation. This gives me a compile error (might have something to do with the elements being of type <code>Eigen::Vector3d</code> or <code>velocity</code> being a class member). Also, I read atomic operations are very <strong>slow</strong> compared to having a private variable for each thread and doing a reduction in the end. So that's what I would like to do, I think.</p>
  
  <p>I have come up with this:</p>
  
  <pre><code>#pragma omp parallel
  {
      // these variables are local to each thread
      std::vector&lt;Eigen::Vector3d&gt; velocity_local(velocity.size());
      std::fill(velocity_local.begin(), velocity_local.end(), Eigen::Vector3d(0,0,0));
  
      #pragma omp for
      for (size_t i = 0; i &lt; clusters.size(); ++i)
      {
          const std::set&lt;int&gt;&amp; cluster = clusters[i];
          // ... expensive calculations ...
          for (int j : cluster)
              velocity_local[j] += f(j); // save results from the previous calculations
      } 
  
      // now each thread can save its results to the global variable
      #pragma omp critical
      {
          for (size_t i = 0; i &lt; velocity_local.size(); ++i)
              velocity[i] += velocity_local[i];
      }
  } 
  </code></pre>
  
  <p>Is this a good solution? Is it the <em>best</em> solution? (Is it even <em>correct</em>?)</p>
  
  <p>Further thoughts: Using the <code>reduce</code> clause (instead of the <code>critical</code> section) throws a compiler error. I think this is because <code>velocity</code> is a class member.</p>
  
  <p>I have tried to find a question with a similar problem, and <a href="https://stackoverflow.com/questions/22073639/openmp-with-parallel-reduction-in-for-loop">this</a> question looks like it's almost the same. But I think my case <em>might</em> differ because the last step includes a <code>for</code> loop. Also the question whether this is the best approach still holds.</p>
  
  <p><strong>Edit:</strong> As request per comment: The <code>reduction</code> clause...</p>
  
  <pre><code>    #pragma omp parallel reduction(+:velocity)
      for (omp_int i = 0; i &lt; velocity_local.size(); ++i)
          velocity[i] += velocity_local[i];
  </code></pre>
  
  <p>...throws the following error:</p>
  
  <p><em>error C3028: 'ShapeMatching::velocity' : only a variable or static data member can be used in a data-sharing clause</em></p>
  
  <p>(similar error with <code>g++</code>)</p>