integrate Hauke's 2x2 direct symmetric eigenvalues solver

Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h

@@ -500,10 +500,10 @@ template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_
 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
 {
   typedef typename SolverType::MatrixType MatrixType;
+  typedef typename SolverType::RealVectorType VectorType;
   typedef typename SolverType::Scalar Scalar;
   
-  template<typename Roots>
+  inline static void computeRoots(const MatrixType& m, VectorType& roots)
-  inline static void computeRoots(const MatrixType& m, Roots& roots)
   {
     using std::sqrt;
     using std::atan2;
@@ -561,7 +561,7 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
     
     MatrixType& eivecs = solver.m_eivec;
-    typename SolverType::RealVectorType& eivals = solver.m_eivalues;
+    VectorType& eivals = solver.m_eivalues;
   
     // map the matrix coefficients to [-1:1] to avoid over- and underflow.
     Scalar scale = mat.cwiseAbs().maxCoeff();
@@ -610,6 +610,91 @@ template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3
   }
 };
 
+// 2x2 direct eigenvalues decomposition, code from Hauke Heibel
+template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false>
+{
+  typedef typename SolverType::MatrixType MatrixType;
+  typedef typename SolverType::RealVectorType VectorType;
+  typedef typename SolverType::Scalar Scalar;
+  
+  inline static void run(SolverType& solver, const MatrixType& mat, int options)
+  {
+    eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
+    eigen_assert((options&~(EigVecMask|GenEigMask))==0
+            && (options&EigVecMask)!=EigVecMask
+            && "invalid option parameter");
+    bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
+    
+    MatrixType& eivecs = solver.m_eivec;
+    VectorType& eivals = solver.m_eivalues;
+  
+    // map the matrix coefficients to [-1:1] to avoid over- and underflow.
+    Scalar scale = mat.cwiseAbs().maxCoeff();
+    scale = std::max(scale,Scalar(1));
+    MatrixType scaledMat = mat / scale;
+
+    // Compute the eigenvalues
+    const Scalar xx = scaledMat(0,0);
+    const Scalar yy = scaledMat(1,1);
+    const Scalar xy = scaledMat(1,0);
+    
+    const Scalar root = std::sqrt((xx-yy)*(xx-yy)+Scalar(4.0)*xy*xy);
+    const Scalar eig_val1 = Scalar(0.5)*(xx+yy+root);
+    const Scalar eig_val2 = Scalar(0.5)*(xx+yy-root);
+    
+    using std::sqrt;
+    using std::abs;
+    
+    // compute the eigen vectors
+    if(computeEigenvectors)
+    {
+      Scalar nx, ny;
+      if (xx-yy > 0)
+      {
+        Scalar len(1.0);
+        if (-yy+xx+root != Scalar(0.0))
+        {
+          ny = -Scalar(2.0)*xy/(-yy+xx+root);
+          len = Scalar(1.0)/sqrt(1+ny*ny);
+          ny *= len;
+        }
+        else
+          ny = Scalar(0.0);
+        nx = -len;
+      }
+      else
+      {
+        Scalar len = Scalar(1.0);
+        if (-yy+xx-root != Scalar(0.0))
+        {       
+          nx = -Scalar(2.0)*xy/(-yy+xx-root);
+          len = Scalar(1.0)/sqrt(1+nx*nx);
+          nx *= len;
+        }
+        else
+          nx = Scalar(0.0);
+        ny = len;
+      }
+      
+      int id0(0), id1(1);
+      if(eig_val1>eig_val2)
+        std::swap(id0,id1);
+
+      eivecs.col(id0) << nx, ny;
+      eivecs.col(id1) << -ny, nx;
+      eivals(id0) = eig_val1;
+      eivals(id1) = eig_val2;
+    }
+    
+    // Rescale back to the original size.
+    eivals *= scale;
+    
+    solver.m_info = Success;
+    solver.m_isInitialized = true;
+    solver.m_eigenvectorsOk = computeEigenvectors;
+  }
+};
+
 }
 
 template<typename MatrixType>