title: Diagonalization of a 2x2 self-adjoined (hermitian) matrix
body: <p>Diagonalizing a 2x2 hermitian matrix is simple, it can be done analytically. However, when it comes to calculating the eigenvalues and eigenvectors over >10^6 times, it is important to do it as efficient as possible. Especially if the off-diagonal elements can vanish it is not possible to use one formula for the eigenvectors: An if-statement is necessary, which of course slows down the code. Thus, I thought using Eigen, where it's stated that the diagonalization of 2x2 and 3x3 matrices is optimized, would be still a good choice:</p>
  
  <p>using</p>
  
  <pre><code>const std::complex&lt;double&gt; I ( 0.,1. );
  inline double block_distr ( double W )
  {
    return (-W/2. + rand() * W/RAND_MAX);
  }
  </code></pre>
  
  <p>a test-loop would be</p>
  
  <pre><code>...
  SelfAdjointEigenSolver&lt;Matrix&lt;complex&lt; double &gt;, 2, 2&gt; &gt; ces;
  Matrix&lt;complex&lt; double &gt;, 2, 2&gt; X;
  
  for (int i = 0 ; i &lt;iter_MAX; ++i) {
    a00=block_distr(100.);
    a11=block_distr(100.);
    re_a01=block_distr(100.);
    im_a01=block_distr(100.);
  
    X(0,0)=a00;
    X(1,0)=re_a01-I*im_a01;
    //only the lower triangular part is referenced! X(0,1)=0.; &lt;--- not necessary
    X(1,1)=a11;
    ces.compute(X,ComputeEigenvectors);
  }
  </code></pre>
  
  <p>Writing the loop without Eigen, using directly the formulas for eigenvalues and eigenvectors of a hermitian matrix and an if-statement to check if the off diagonal is zero, is a factor of 5 faster. Am I not using Eigen properly or is such an overhead normal? Are there other lib.s which are optimized for small self-adjoint matrices?</p>