title: C++ - Fastor vs. Eigen vs C-tran: Speed differences in tensor contraction
body: <p>I am currently looking for a C++ library (preferably header only) for higher order tensor contractions with einstein summation like notation.</p>
  
  <p>Fastor (<a href="https://github.com/romeric/Fastor" rel="nofollow noreferrer">https://github.com/romeric/Fastor</a>) seems to be perfectly suited, and since Eigen (which I am using a lot) has a tensor module, I made a small comparison, including a benchmark of a simple loop implementation:</p>
  
  <pre><code>#include &lt;Fastor.h&gt;
  #include "unsupported/Eigen/CXX11/Tensor"
  #include &lt;ctime&gt;
  
  int main() {
  clock_t tstart, tend; 
   {
      using namespace Fastor;
      Tensor&lt;double,20,50,10&gt; A; 
      Tensor&lt;double,50,10,20,4&gt; B;
      Tensor&lt;double, 4&gt; C;
      Tensor&lt;double, 10, 10, 4&gt; D;
  
      A.ones();
      B.ones();
      C.zeros();
      D.zeros();
      enum {I,J,K,L,M,N};
  
      tstart = clock();
  
      C += einsum&lt;Index&lt;I,J,K&gt;,Index&lt;J,K,I,M&gt;&gt;(A,B);
  
      tend = clock(); 
      std::cout &lt;&lt; "Fastor, C_M = A_IJK * B_JKIM:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
  
      tstart = clock();
  
      D += einsum&lt;Index&lt;I,J,K&gt;,Index&lt;J,M,I,N&gt;&gt;(A,B);
  
      tend = clock(); 
      std::cout &lt;&lt; "Fastor, D_KMN = A_IJ * B_JMIN:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
   }
  
   {
       using namespace Eigen;
  
      TensorFixedSize&lt;double, Sizes&lt;20, 50, 10&gt;&gt; A;
      TensorFixedSize&lt;double, Sizes&lt;50,10, 20, 4&gt;&gt; B;
      TensorFixedSize&lt;double, Sizes&lt;4&gt;&gt; C;
      TensorFixedSize&lt;double, Sizes&lt;10,10,4&gt;&gt; D;
  
      A.setConstant(1);
      B.setConstant(1);
      C.setConstant(0);
      D.setConstant(0);
  
       array&lt;IndexPair&lt;int&gt;,3&gt; IJK_JKIM = {
           IndexPair&lt;int&gt;(0, 2),
           IndexPair&lt;int&gt;(1, 0),
           IndexPair&lt;int&gt;(2, 1),
       };
  
       array&lt;IndexPair&lt;int&gt;,2&gt; IJK_JMIN = {
           IndexPair&lt;int&gt;(0, 2),
           IndexPair&lt;int&gt;(1, 0),
       };
  
      tstart = clock();
      C += A.contract(  B,  IJK_JKIM) ;
      tend = clock(); 
  
      std::cout &lt;&lt; "Eigen, C_M = A_IJK * B_JKIM:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
  
      tstart = clock();
       D += A.contract(  B,  IJK_JMIN) ;
      tend = clock(); 
  
      std::cout &lt;&lt; "Eigen, D_KMN = A_IJ * B_JMIN:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
  
  
   }
  
   {
       using namespace Eigen;
  
       double A [20][50][10]  = {1} ;
       double B [50][10][20][4]  = {1};
       double C [4]  = {};
       double D [10][10][4]  = {};
  
      for ( int i = 0; i &lt; 20; i++)
          for ( int j = 0; j &lt; 50; j++)
              for ( int k = 0; k &lt; 10; k++)
                  A[i][j][k] = 1.0;
  
      for ( int i = 0; i &lt; 20; i++)
          for ( int j = 0; j &lt; 50; j++)
              for ( int k = 0; k &lt; 10; k++)
                  for ( int l = 0; l &lt; 4; l++)
                      B[j][k][i][l] = 1.0;
  
  
      tstart = clock();
  
      for ( int i = 0; i &lt; 20; i++)
          for ( int j = 0; j &lt; 50; j++)
              for ( int k = 0; k &lt; 10; k++)
                  for ( int l = 0; l &lt; 4; l++)
                      C[l] += A[i][j][k] * B [j][k][i][l];
  
      tend = clock(); 
      std::cout &lt;&lt; "CTran, C_M = A_IJK * B_JKIM:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
  
      tstart = clock();
      for ( int i = 0; i &lt; 20; i++)
          for ( int j = 0; j &lt; 50; j++)
              for ( int k = 0; k &lt; 10; k++)
                  for ( int m = 0; m &lt; 10; m++)
                      for ( int n = 0; n &lt; 4; n++)
                          D[k][m][n] += A[i][j][k] * B [j][m][i][n];
  
      tend = clock(); 
  
      std::cout &lt;&lt; "CTran, D_KMN = A_IJ * B_JMIN:\t"&lt;&lt; tend - tstart &lt;&lt; std::endl;
  
   }
  
  return 0;
  }
  </code></pre>
  
  <p>Output for me is:</p>
  
  <pre><code>Fastor, C_M = A_IJK * B_JKIM:   206
  Fastor, D_KMN = A_IJ * B_JMIN:  2230
  Eigen, C_M = A_IJK * B_JKIM:    1286
  Eigen, D_KMN = A_IJ * B_JMIN:   898
  CTran, C_M = A_IJK * B_JKIM:    2
  CTran, D_KMN = A_IJ * B_JMIN:   2
  </code></pre>
  
  <p>Compiled with g++ 9.1.0 (Arch Linux) using</p>
  
  <pre><code>g++ test.cpp -O3 -std=c++17  -I Fastor -isystem eigen-eigen-323c052e1731 -o test
  </code></pre>
  
  <p>So it seems that Fastor is considerably faster for example 1 than Eigen, but slower for example 2.
  However, both libraries are incredibely slower than the naive loop implementation!
  Is there an explanation for these contradicting results? Thank you in advance!</p>