title: Why are elementwise additions much faster in separate loops than in a combined loop?
body: <p>Suppose <code>a1</code>, <code>b1</code>, <code>c1</code>, and <code>d1</code> point to heap memory, and my numerical code has the following core loop.</p>
  <pre><code>const int n = 100000;
  
  for (int j = 0; j &lt; n; j++) {
      a1[j] += b1[j];
      c1[j] += d1[j];
  }
  </code></pre>
  <p>This loop is executed 10,000 times via another outer <code>for</code> loop. To speed it up, I changed the code to:</p>
  <pre><code>for (int j = 0; j &lt; n; j++) {
      a1[j] += b1[j];
  }
  
  for (int j = 0; j &lt; n; j++) {
      c1[j] += d1[j];
  }
  </code></pre>
  <p>Compiled on <a href="http://en.wikipedia.org/wiki/Visual_C++#32-bit_versions" rel="noreferrer">Microsoft Visual C++ 10.0</a> with full optimization and <a href="http://en.wikipedia.org/wiki/SSE2" rel="noreferrer">SSE2</a> enabled for 32-bit on a <a href="http://en.wikipedia.org/wiki/Intel_Core_2" rel="noreferrer">Intel Core 2</a> Duo (x64), the first example takes 5.5 seconds and the double-loop example takes only 1.9 seconds.</p>
  <p>Disassembly for the first loop basically looks like this (this block is repeated about five times in the full program):</p>
  <pre><code>movsd       xmm0,mmword ptr [edx+18h]
  addsd       xmm0,mmword ptr [ecx+20h]
  movsd       mmword ptr [ecx+20h],xmm0
  movsd       xmm0,mmword ptr [esi+10h]
  addsd       xmm0,mmword ptr [eax+30h]
  movsd       mmword ptr [eax+30h],xmm0
  movsd       xmm0,mmword ptr [edx+20h]
  addsd       xmm0,mmword ptr [ecx+28h]
  movsd       mmword ptr [ecx+28h],xmm0
  movsd       xmm0,mmword ptr [esi+18h]
  addsd       xmm0,mmword ptr [eax+38h]
  </code></pre>
  <p>Each loop of the double loop example produces this code (the following block is repeated about three times):</p>
  <pre><code>addsd       xmm0,mmword ptr [eax+28h]
  movsd       mmword ptr [eax+28h],xmm0
  movsd       xmm0,mmword ptr [ecx+20h]
  addsd       xmm0,mmword ptr [eax+30h]
  movsd       mmword ptr [eax+30h],xmm0
  movsd       xmm0,mmword ptr [ecx+28h]
  addsd       xmm0,mmword ptr [eax+38h]
  movsd       mmword ptr [eax+38h],xmm0
  movsd       xmm0,mmword ptr [ecx+30h]
  addsd       xmm0,mmword ptr [eax+40h]
  movsd       mmword ptr [eax+40h],xmm0
  </code></pre>
  <p>The question turned out to be of no relevance, as the behavior severely depends on the sizes of the arrays (n) and the CPU cache. So if there is further interest, I rephrase the question:</p>
  <ul>
  <li><p>Could you provide some solid insight into the details that lead to the different cache behaviors as illustrated by the five regions on the following graph?</p>
  </li>
  <li><p>It might also be interesting to point out the differences between CPU/cache architectures, by providing a similar graph for these CPUs.</p>
  </li>
  </ul>
  <p>Here is the full code. It uses <a href="https://www.threadingbuildingblocks.org/" rel="noreferrer">TBB</a> <code>Tick_Count</code> for higher resolution timing, which can be disabled by not defining the <code>TBB_TIMING</code> Macro:</p>
  <pre><code>#include &lt;iostream&gt;
  #include &lt;iomanip&gt;
  #include &lt;cmath&gt;
  #include &lt;string&gt;
  
  //#define TBB_TIMING
  
  #ifdef TBB_TIMING   
  #include &lt;tbb/tick_count.h&gt;
  using tbb::tick_count;
  #else
  #include &lt;time.h&gt;
  #endif
  
  using namespace std;
  
  //#define preallocate_memory new_cont
  
  enum { new_cont, new_sep };
  
  double *a1, *b1, *c1, *d1;
  
  
  void allo(int cont, int n)
  {
      switch(cont) {
        case new_cont:
          a1 = new double[n*4];
          b1 = a1 + n;
          c1 = b1 + n;
          d1 = c1 + n;
          break;
        case new_sep:
          a1 = new double[n];
          b1 = new double[n];
          c1 = new double[n];
          d1 = new double[n];
          break;
      }
  
      for (int i = 0; i &lt; n; i++) {
          a1[i] = 1.0;
          d1[i] = 1.0;
          c1[i] = 1.0;
          b1[i] = 1.0;
      }
  }
  
  void ff(int cont)
  {
      switch(cont){
        case new_sep:
          delete[] b1;
          delete[] c1;
          delete[] d1;
        case new_cont:
          delete[] a1;
      }
  }
  
  double plain(int n, int m, int cont, int loops)
  {
  #ifndef preallocate_memory
      allo(cont,n);
  #endif
  
  #ifdef TBB_TIMING   
      tick_count t0 = tick_count::now();
  #else
      clock_t start = clock();
  #endif
          
      if (loops == 1) {
          for (int i = 0; i &lt; m; i++) {
              for (int j = 0; j &lt; n; j++){
                  a1[j] += b1[j];
                  c1[j] += d1[j];
              }
          }
      } else {
          for (int i = 0; i &lt; m; i++) {
              for (int j = 0; j &lt; n; j++) {
                  a1[j] += b1[j];
              }
              for (int j = 0; j &lt; n; j++) {
                  c1[j] += d1[j];
              }
          }
      }
      double ret;
  
  #ifdef TBB_TIMING   
      tick_count t1 = tick_count::now();
      ret = 2.0*double(n)*double(m)/(t1-t0).seconds();
  #else
      clock_t end = clock();
      ret = 2.0*double(n)*double(m)/(double)(end - start) *double(CLOCKS_PER_SEC);
  #endif
      
  #ifndef preallocate_memory
      ff(cont);
  #endif
  
      return ret;
  }
  
  
  void main()
  {   
      freopen(&quot;C:\\test.csv&quot;, &quot;w&quot;, stdout);
  
      char *s = &quot; &quot;;
  
      string na[2] ={&quot;new_cont&quot;, &quot;new_sep&quot;};
  
      cout &lt;&lt; &quot;n&quot;;
  
      for (int j = 0; j &lt; 2; j++)
          for (int i = 1; i &lt;= 2; i++)
  #ifdef preallocate_memory
              cout &lt;&lt; s &lt;&lt; i &lt;&lt; &quot;_loops_&quot; &lt;&lt; na[preallocate_memory];
  #else
              cout &lt;&lt; s &lt;&lt; i &lt;&lt; &quot;_loops_&quot; &lt;&lt; na[j];
  #endif
              
      cout &lt;&lt; endl;
  
      long long nmax = 1000000;
  
  #ifdef preallocate_memory
      allo(preallocate_memory, nmax);
  #endif
      
      for (long long n = 1L; n &lt; nmax; n = max(n+1, long long(n*1.2)))
      {
          const long long m = 10000000/n;
          cout &lt;&lt; n;
  
          for (int j = 0; j &lt; 2; j++)
              for (int i = 1; i &lt;= 2; i++)
                  cout &lt;&lt; s &lt;&lt; plain(n, m, j, i);
          cout &lt;&lt; endl;
      }
  }
  </code></pre>
  <p>It shows FLOP/s for different values of <code>n</code>.</p>
  <p><a href="https://i.stack.imgur.com/426bP.gif" rel="noreferrer"><img src="https://i.stack.imgur.com/426bP.gif" alt="Performace chart" /></a></p>