title: How does #pragma simd reduction(<operator>:<variable>) work under the hood?
body: <p>I would like to know in more detail how the simd reduction clause used by Intel compilers works under the hood.</p>
  
  <p>In particular, for a loop of the form</p>
  
  <pre><code>double x = x_initial;
  #pragma simd reduction(&lt;operator1&gt;:x)
  for( int i = 0; i &lt; N; i++ )
    x &lt;operator2&gt; some_value;
  </code></pre>
  
  <p>my naive guess is as follows: 
  The compiler initializes a private copy of x for each vector lane, then iterates through the loop one vector width at a time.  If the vector width is 4 doubles, for example, this would correspond to N/4 iterations plus a peel loop at the end.  At each step of the iteration, each lane's private copy of x is updated using <code>operator2</code>, then at the end, the 4 vector lanes' private copies are combined using <code>operator1</code>.  The <a href="https://software.intel.com/sites/default/files/8c/a9/CompilerAutovectorizationGuide.pdf" rel="noreferrer" title="auto-vectorization guide">auto-vectorization guide</a> does not appear to address this directly.</p>
  
  <p>I did some experimentation and found some results that agree with my expectation and some that don't.  For example, I tried the case</p>
  
  <pre><code>double x = 1;
  #pragma simd reduction(*:x) assert
  for( int i = 0; i &lt; 16; i++ )
    x += a[i];  // All elements of a are equal to 3.0
  cout &lt;&lt; "x after (*:x), x += a[i] loop:  " &lt;&lt; x &lt;&lt; endl; 
  </code></pre>
  
  <p>where <code>operator1</code> is * and <code>operator2</code> is +=.  When I compile for avx2, which has a vector width of 4 doubles, the output is 28561 = ( 1 + 4*a[i] )^4.  This implies that the code first initializes 4 lane-private copies of x to 1, then adds 3 to each of those copies 4 times as the 4-double-wide vector lane iterates across the trip count of 16.  Each lane-private copy of x is now equal to 13.  Finally, the lane-private copies are combined (reduced) using <code>operator2</code> which is *, yielding 13*13*13*13 = 28561.</p>
  
  <p>However, when I switch the * and + operators, like so</p>
  
  <pre><code>x = 1;
  #pragma simd reduction(+:x) assert
  for( int i = 0; i &lt; 16; i++ )
    x *= a[i];
  cout &lt;&lt; "x after (+:x), x *= a[i] loop:  " &lt;&lt; x &lt;&lt; endl;
  </code></pre>
  
  <p>and compile again for avx2, the output is 1.0.  If my theory were correct, each vector lane should end up containing a value of 1*3^4, which would then be combined using + to yield 4*3^4 = 324.  Evidently this is not the case.  What am I missing?</p>