title: Eigen C++ running slower when compiled with MKL
body: <p>I have recently started to use Eigen (version 3.3.1), running a benchmark against Armadillo on a simple matrix operation at the core of an OLS regression, that is computing the inverse of the product of a matrix by itself, I noticed that Eigen was running slower when compiled with MKL library than without it for this kind of operation. I was wondering if my compilation instructions were wrong. I also tried to realize this operation calling directly MKL BLAS and LAPACK routines and got a much faster result, as fast as Armadillo. I cannot explain such a poor performance especially for float type.</p>
  
  <p>I wrote the code below to realise this benchmark:</p>
  
  <pre><code>#define ARMA_DONT_USE_WRAPPER
  #define ARMA_NO_DEBUG
  #include &lt;armadillo&gt;
  
  #define EIGEN_NO_DEBUG
  #define EIGEN_NO_STATIC_ASSERT
  #define EIGEN_USE_MKL_ALL
  #include &lt;Eigen/Dense&gt;
  
  template &lt;typename T&gt;
  using Matrix = Eigen::Matrix&lt;T, Eigen::Dynamic, Eigen::Dynamic&gt;;
  
  #ifdef USE_FLOAT
  using T = float;
  #else
  using T = double;
  #endif
  
  int main()
  {
      arma::wall_clock timer;
  
      int niter = 1000000;
      int n = 1000;
      int k = 20;
  
      arma::Mat&lt;T&gt; Xa = arma::cumsum(arma::randn&lt;arma::Mat&lt;T&gt;&gt;(n, k));
      Matrix&lt;T&gt; Xe = Matrix&lt;T&gt;::Map(Xa.memptr(), Xa.n_rows, Xa.n_cols);
  
      // Armadillo compiled with MKL
      timer.tic();
      for (int i = 0; i &lt; niter; ++i) {
          arma::Mat&lt;T&gt; iX2a = (Xa.t() * Xa).i();
      }
      std::cout &lt;&lt; "...Elapsed time: " &lt;&lt; timer.toc() &lt;&lt; "\n";
  
      // Eigen compiled with MKL
      timer.tic();
      for (int i = 0; i &lt; niter; ++i) {
          Matrix&lt;T&gt; iX2e = (Xe.transpose() * Xe).inverse();
      }
      std::cout &lt;&lt; "...Elapsed time: " &lt;&lt; timer.toc() &lt;&lt; "\n";*/
  
      // Eigen Matrix with MKL routines
      timer.tic();
      for (int i = 0; i &lt; niter; ++i) {
          Matrix&lt;T&gt; iX2e =  Matrix&lt;T&gt;::Zero(k, k);
          // first stage =&gt; computing square matrix trans(X) * X
          #ifdef USE_FLOAT
          cblas_ssyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &amp;Xe(0,0), n, 0.0, &amp;iX2e(0,0), k);
          #else
          cblas_dsyrk(CblasColMajor, CblasLower, CblasTrans, k, n, 1.0, &amp;Xe(0,0), n, 0.0, &amp;iX2e(0,0), k);
          #endif
          // getting upper part  
          for (int i = 0; i &lt; k; ++i)
              for (int j = i + 1; j &lt; k; ++j)
                  iX2e(i, j) = iX2e(j, i);
          // second stage =&gt; inverting square matrix
          // initializing pivots
          int* ipiv = new int[k];
          // factorizing matrix
          #ifdef USE_FLOAT 
          LAPACKE_sgetrf(LAPACK_COL_MAJOR, k, k, &amp;iX2e(0,0), k, ipiv);
          #else
          LAPACKE_dgetrf(LAPACK_COL_MAJOR, k, k, &amp;iX2e(0,0), k, ipiv); 
          #endif
          // computing the matrix inverse
          #ifdef USE_FLOAT 
          LAPACKE_sgetri(LAPACK_COL_MAJOR, k, &amp;iX2e(0,0), k, ipiv);
          #else
          LAPACKE_dgetri(LAPACK_COL_MAJOR, k, &amp;iX2e(0,0), k, ipiv);
          #endif
          delete[] ipiv;
      }
      std::cout &lt;&lt; "...Elapsed time: " &lt;&lt; timer.toc() &lt;&lt; "\n";
  }
  </code></pre>
  
  <p>I compile this file called test.cpp with:</p>
  
  <p>g++ -std=c++14 -Wall -O3 -march=native -DUSE_FLOAT test.cpp -o run -L${MKLROOT}/lib/intel64 -Wl,--no-as-needed -lmkl_gf_lp64 -lmkl_sequential -lmkl_core</p>
  
  <p>I get the following results (on Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz)</p>
  
  <ul>
  <li>for double type:</li>
  </ul>
  
  <p>Armadillo with MKL => 64.0s</p>
  
  <p>Eigen with MKL => 72.2s</p>
  
  <p>Eigen alone => 68.7s</p>
  
  <p>Pure MKL => 64.9s  </p>
  
  <ul>
  <li>for float type:</li>
  </ul>
  
  <p>Armadillo with MKL => 38.2s</p>
  
  <p>Eigen with MKL => 61.1s</p>
  
  <p>Eigen alone => 42.6s</p>
  
  <p>Pure MKL => 38.9s</p>
  
  <p>NB: I run this test for a project which will not use very large matrix, I do not need parallelization at this level, my biggest matrix will probably be 2000 rows for 25 columns, moreover I will need to go parallel at a higher level so I want to avoid any kind of nested parallelism which could slow down my code.</p>