title: Understanding SSE3 Matrix Multiplication Optimization
body: <p>With reference to <a href="http://blogs.msdn.com/b/xiangfan/archive/2009/04/28/optimize-your-code-matrix-multiplication.aspx" rel="nofollow">http://blogs.msdn.com/b/xiangfan/archive/2009/04/28/optimize-your-code-matrix-multiplication.aspx</a>. </p>
  
  <pre><code>template&lt;&gt;
  void SeqMatrixMult4(int size, float** m1, float** m2, float** result)
  {
      Transpose(size, m2);
      for (int i = 0; i &lt; size; i++) {
          for (int j = 0; j &lt; size; j++) {
              __m128 c = _mm_setzero_ps();
  
              for (int k = 0; k &lt; size; k += 4) {
                  c = _mm_add_ps(c, _mm_mul_ps(_mm_load_ps(&amp;m1[i][k]), _mm_load_ps(&amp;m2[j][k])));
              }
              c = _mm_hadd_ps(c, c);
              c = _mm_hadd_ps(c, c);
              _mm_store_ss(&amp;result[i][j], c);
          }
      }
      Transpose(size, m2);
  }
  </code></pre>
  
  <p>Why is there 2 more <code>_mm_hadd_ps(c, c)</code> after the inner most for loop? To verify my understanding: this code loads 4 floats from m1 and another 4 from m2, then multiplies them resulting in 4 floats (<code>__m128</code>). Then I sum them into <code>c</code> (at this point, its still 4 floats?). Then after the for loop I <code>hadd</code> this result twice? What does that do?  </p>
  
  <hr>
  
  <p>My code slightly re-written produces the wrong result it appears</p>
  
  <pre><code>long long start, end;
  __m128 v1, v2, vMul, vRes;
  vRes = _mm_setzero_ps();
  
  start = wall_clock_time();
  transpose_matrix(m2);
  for (int i = 0; i &lt; SIZE; i++) {
      for (int j = 0; j &lt; SIZE; j++) {
          float tmp = 0;
          for (int k = 0; k &lt; SIZE; k+=4) {
              v1 = _mm_load_ps(&amp;m1[i][k]);
              v2 = _mm_load_ps(&amp;m2[j][k]);
              vMul = _mm_mul_ps(v1, v2);
  
              vRes = _mm_add_ps(vRes, vMul);
          }
          vRes = _mm_hadd_ps(vRes, vRes);
          _mm_store_ss(&amp;result[i][j], vRes);
      }
  }
  end = wall_clock_time();
  fprintf(stderr, "Optimized Matrix multiplication took %1.2f seconds\n", ((float)(end - start))/1000000000);
  
  // reverse the transposition
  transpose_matrix(m2);
  </code></pre>