title: Matrix traversal optimization
body: <p>Given a <code>n x n</code> matrix of <code>ints</code>, I have an algorithm that at each step of a for loop of range <code>n</code> traverses and modifies the matrix. Here is the code:</p>
  <pre class="lang-c prettyprint-override"><code>typedef int **Matrix;
  
  void floyd_slow(Matrix dist, int n)
  {
      int d;
  
      for (int k=0; k&lt;n; k++)
      {
          for (int i=0; i&lt;n; i++)
          {
              for (int j=0; j&lt;n; j++)
                  if ((d=dist[k][j]+dist[i][k])&lt;dist[i][j])
                      dist[i][j]=d;
          }
      }
      for (int i=0; i&lt;n; i++)
          dist[i][i]=0;
  }
  </code></pre>
  <p>The matrix is built as an array of <code>n*n</code> ints and for each index line <code>i</code>, <code>dist[i]</code> is the address of the row of index <code>i</code> [the above code is the <a href="https://rosettacode.org/wiki/Floyd-Warshall_algorithm#C" rel="nofollow noreferrer">standard way</a> to write the Floyd-Warshall algorithm but my question is not about this algorithm by itself].</p>
  <p>The following drawing tries to explain how the matrix is processed:</p>
  <p><a href="https://i.stack.imgur.com/3IruW.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3IruW.png" alt="matrix processing" /></a></p>
  <p>At each step of the loop of index <code>k</code>, the underlying matrix is traversed line by line.</p>
  <p>Now, consider the following transformation of the previous code:</p>
  <pre class="lang-c prettyprint-override"><code>void relax(Matrix dist, int n, int* rowk, int* colk)
  {
      int d;
      for (int i=0; i&lt;n; i++)
          for (int j=0; j&lt;n; j++)
              if ((d=rowk[j]+colk[i])&lt;dist[i][j])
                  dist[i][j]=d;
  }
  
  void floyd_fast(Matrix dist, int n)
  {
      int i, k;
  
      int* colk=malloc(n*sizeof(int));
  
      if (!colk)
          exit(EXIT_FAILURE);
  
      for (k=0; k&lt;n; k++)
      {
          int* rowk =dist[k];
  
          for (i=0; i&lt;n; i++)
              colk[i]=dist[i][k];
  
          relax(dist, n, rowk, colk);
      }
      free(colk);
      for (i=0; i&lt;n; i++)
          dist[i][i]=0;
  }
  </code></pre>
  <p>At every step, the elements of the matrix are accessed in the same order as in the previous algorithm.</p>
  <p>The only difference is that at each step <code>k</code> in the exterior loop, the column of index <code>k</code> is <strong>copied</strong> into a temporary array, cf. the <code>colk</code> malloc above. It results that the element at position <code>(i, k)</code> is read from this array instead of being accessed directly from the matrix.</p>
  <p>This innocuous change leads in fact to a significant speedup: you gain a factor 4 if n=1000.</p>
  <p>I know that in C, it's faster to traverse an array in row major order but this is always the case here. So i was wondering why there is a speedup so important. Is it related to cache optimisation?</p>
  <h1>Complete code</h1>
  <pre class="lang-c prettyprint-override"><code>#include &lt;stdio.h&gt;
  #include &lt;stdlib.h&gt;
  #include &lt;time.h&gt;
  #include &lt;string.h&gt;
  
  typedef int **Matrix;
  
  void floyd_slow(Matrix dist, int n)
  {
      int d;
  
      for (int k=0; k&lt;n; k++)
      {
          for (int i=0; i&lt;n; i++)
          {
              for (int j=0; j&lt;n; j++)
                  if ((d=dist[k][j]+dist[i][k])&lt;dist[i][j])
                      dist[i][j]=d;
          }
      }
      for (int i=0; i&lt;n; i++)
          dist[i][i]=0;
  }
  
  void relax(Matrix dist, int n, int* rowk, int* colk)
  {
      int d;
      for (int i=0; i&lt;n; i++)
          for (int j=0; j&lt;n; j++)
              if ((d=rowk[j]+colk[i])&lt;dist[i][j])
                  dist[i][j]=d;
  }
  
  void floyd_fast(Matrix dist, int n)
  {
      int i, k;
  
      int* colk=malloc(n*sizeof(int));
  
      if (!colk)
          exit(EXIT_FAILURE);
  
      for (k=0; k&lt;n; k++)
      {
          int* rowk =dist[k];
  
          for (i=0; i&lt;n; i++)
              colk[i]=dist[i][k];
  
          relax(dist, n, rowk, colk);
      }
      free(colk);
      for (i=0; i&lt;n; i++)
          dist[i][i]=0;
  }
  
  
  void print(Matrix dist, int n)
  {
      int i, j;
  
      for (i=0; i&lt;n; i++)
      {
          for (j=0; j&lt;n; j++)
              printf(&quot;%d &quot;, dist[i][j]);
          printf(&quot;\n&quot;);
      }
  }
  
  void test_slow(Matrix dist, int n)
  {
      clock_t now=clock();
      floyd_slow(dist, n);
  //         print(dist, n);
      int *p=dist[0];
      free(dist);
      free(p);
  
      fprintf(stderr, &quot;Elapsed slow: %.2f s\n&quot;,
              (double) (clock() - now) / CLOCKS_PER_SEC);
  }
  
  void test_fast(Matrix dist, int n)
  {
      clock_t now=clock();
      floyd_fast(dist, n);
  //         print(dist, n);
      int *p=dist[0];
      free(dist);
      free(p);
      
      fprintf(stderr, &quot;Elapsed fast: %.2f s\n&quot;,
              (double) (clock() - now) / CLOCKS_PER_SEC);
  }
  
  int * data(int n)
  {
      int N=n*n;
      int *t=malloc(N*sizeof(int));
      if (!t)
          exit(EXIT_FAILURE);  
  
      srand(time(NULL)); 
      
      for (int i=0; i&lt;N;i++)
          t[i]=(1+rand())%10;
          
      return t;    
  }
  
  Matrix getMatrix(int *t, int n)
  {
      int N=n*n;
      int *tt=malloc(N*sizeof(int));
      Matrix mat=malloc(n*sizeof(int*));
      
      if (!tt || !mat)
          exit(EXIT_FAILURE); 
  
      memcpy(tt, t, N*sizeof(int));
      
      for (int i=0; i&lt;n;i++)
          mat[i]=&amp;tt[i*n];
      
      return mat;
  }
  
  
  int main(void)
  {
      int n=1000;
      int *t=data(n);
      
      Matrix mat_slow=getMatrix(data(n), n);
      Matrix mat_fast=getMatrix(data(n), n); 
      
      test_slow(mat_slow, n);
      test_fast(mat_fast, n);
  
      return 0;
  }
  </code></pre>
  <p>Output:</p>
  <pre><code>Elapsed slow: 0.58 s
  Elapsed fast: 0.14 s
  </code></pre>
  <p>Compilation options:</p>
  <pre class="lang-sh prettyprint-override"><code>rm floyd
  gcc -Wall -O3 -march=native -ffast-math -Wno-unused-result -Wno-unused-variable -Wno-unused-but-set-variable -Wno-unused-parameter floyd.c -o floyd -lm
  </code></pre>