title: how to optimize matrix multiplication (matmul) code to run fast on a single processor core
body: <p>I am working on parallel programming concepts and trying to optimize matrix multiplication example on single core. The fastest implementation I came up so far is the following:</p>
  
  <pre><code>/* This routine performs a dgemm operation
   *  C := C + A * B
   * where A, B, and C are lda-by-lda matrices stored in column-major format.
   * On exit, A and B maintain their input values. */    
  void square_dgemm (int n, double* A, double* B, double* C)
  {
    /* For each row i of A */
    for (int i = 0; i &lt; n; ++i)
      /* For each column j of B */
      for (int j = 0; j &lt; n; ++j) 
      {
        /* Compute C(i,j) */
        double cij = C[i+j*n];
        for( int k = 0; k &lt; n; k++ )
            cij += A[i+k*n] * B[k+j*n];
        C[i+j*n] = cij;
      }
  }
  </code></pre>
  
  <p>The results are like below. how to reduce the loops and increase the performance </p>
  
  <pre><code>login4.stampede(72)$ tail -f job-naive.stdout
  Size: 480       Mflop/s:  1818.89       Percentage: 18.95
  Size: 511       Mflop/s:  2291.73       Percentage: 23.87
  Size: 512       Mflop/s:  937.061       Percentage:  9.76
  Size: 639       Mflop/s:  293.434       Percentage:  3.06
  Size: 640       Mflop/s:  270.238       Percentage:  2.81
  Size: 767       Mflop/s:  240.209       Percentage:  2.50
  Size: 768       Mflop/s:  242.118       Percentage:  2.52
  Size: 769       Mflop/s:  240.173       Percentage:  2.50
  Average percentage of Peak = 22.0802
  Grade = 33.1204
  </code></pre>