title: Concatenating Eigen Tensors
body: <p>I am working with tensors in Eigen and need to join two tensors together to form a higher order tensor. As an example:</p>
  
  <p>C_ijkl = A_ik B_jl</p>
  
  <p>I know I can do this with some for loops but is there a better way?</p>
  
  <p>Edit:
  Following the answer below I implemented this into code. Using the g++ compiler with -O3 and comparing an explicit loop versus the contraction shown below I got the following results for 10,000 iterations:</p>
  
  <pre><code>Average time (ns) for explicit loop: 196.782
  Average time (ns) for contraction  : 790.439
  </code></pre>
  
  <p>So it seems like with my test setup a naive explicit loop is faster (though messier) than a contraction.</p>
  
  <p>The code:</p>
  
  <pre><code>#include &lt;iostream&gt;
  #include &lt;Eigen/Core&gt;
  #include &lt;Eigen/Dense&gt;
  #include &lt;unsupported/Eigen/CXX11/Tensor&gt;
  
  #include&lt;chrono&gt;
  typedef std::chrono::high_resolution_clock Clock;
  
  using namespace Eigen;
  
  int main() {
  
  //Initialize
  Tensor&lt;double, 2&gt; A_ij(4, 4);
  Tensor&lt;double, 2&gt; B_kl(4, 4);
  Tensor&lt;double, 4&gt; C_ijkl(4, 4, 4, 4);
  
  A_ij.setValues({{ 1, 2, 3, 4},
                  { 5, 6, 7, 8},
                  { 9,10,11,12},
                  {13,14,15,16}});
  
  B_kl.setValues({{17,18,19,20},
                  {21,22,23,24},
                  {25,26,27,28},
                  {29,30,31,32}});
  
  Tensor&lt;double, 4&gt; Ctrue_ijkl(4, 4, 4, 4);
  
  int nrepeat = 10000;
  double avg_time_explicit;
  double avg_time_contract;
  
  auto t1 = Clock::now();
  auto t2 = Clock::now();
  
  array&lt;Eigen::IndexPair&lt;long&gt;,0&gt; empty_index_list = {};
  
  for (int r=0; r&lt;nrepeat; r++){
      t1 = Clock::now();
      for (int i=0; i&lt;4; i++){
          for (int j=0; j&lt;4; j++){
              for (int k=0; k&lt;4; k++){
                  for (int l=0; l&lt;4; l++){
                      Ctrue_ijkl(i,j,k,l) = A_ij(i,j)*B_kl(k,l);
                  }
              }
          }
      }
      t2 = Clock::now();
      avg_time_explicit += std::chrono::duration_cast&lt;std::chrono::nanoseconds&gt;(t2-t1).count()/static_cast&lt;double&gt;(nrepeat);
  }
  
  
  
  for (int r=0; r&lt;nrepeat; r++){
      t1 = Clock::now();
      C_ijkl = A_ij.contract(B_kl, empty_index_list);
      t2 = Clock::now();
      avg_time_contract += std::chrono::duration_cast&lt;std::chrono::nanoseconds&gt;(t2-t1).count()/static_cast&lt;double&gt;(nrepeat);
  }
  
  //Check to make sure the two methods produce the same result
  Map&lt;VectorXd&gt; mt(C_ijkl.data(), C_ijkl.size());
  Map&lt;VectorXd&gt; mr(Ctrue_ijkl.data(), Ctrue_ijkl.size());
  
  std::cout &lt;&lt; "Ctrue_ijkl == C_ijkl: " &lt;&lt; mt.isApprox(mr) &lt;&lt; "\n";
  
  std::cout &lt;&lt; "Average time (ns) for explicit loop: " &lt;&lt; avg_time_explicit &lt;&lt; "\n";
  std::cout &lt;&lt; "Average time (ns) for contraction:   " &lt;&lt; avg_time_contract &lt;&lt; "\n";
  
  }
  </code></pre>