From cd8b996f99de67035a0504cbaf0a627fb68f0f1d Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Mon, 8 Jun 2015 16:16:42 +0200
Subject: [PATCH] Extend unit test and documentation of
 SelfAdjointEigenSolver::computeDirect

---
 .../src/Eigenvalues/SelfAdjointEigenSolver.h  | 10 ++++--
 test/eigensolver_selfadjoint.cpp              | 32 +++++++++++++------
 2 files changed, 29 insertions(+), 13 deletions(-)

diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 1830b99c1..27a014a96 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -198,17 +198,21 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     EIGEN_DEVICE_FUNC
     SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors);
     
-    /** \brief Computes eigendecomposition of given matrix using a direct algorithm
+    /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm
       *
       * This is a variant of compute(const MatrixType&, int options) which
       * directly solves the underlying polynomial equation.
       * 
-      * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d).
+      * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d).
       * 
-      * This method is usually significantly faster than the QR algorithm
+      * This method is usually significantly faster than the QR iterative algorithm
       * but it might also be less accurate. It is also worth noting that
       * for 3x3 matrices it involves trigonometric operations which are
       * not necessarily available for all scalar types.
+      * 
+      * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues:
+      *   - double: 1e-8
+      *   - float:  1e-3
       *
       * \sa compute(const MatrixType&, int options)
       */
diff --git a/test/eigensolver_selfadjoint.cpp b/test/eigensolver_selfadjoint.cpp
index 6750c7609..41b6d99ab 100644
--- a/test/eigensolver_selfadjoint.cpp
+++ b/test/eigensolver_selfadjoint.cpp
@@ -18,22 +18,34 @@ template<typename MatrixType> void selfadjointeigensolver_essential_check(const
 {
   typedef typename MatrixType::Scalar Scalar;
   typedef typename NumTraits<Scalar>::Real RealScalar;
-  RealScalar largerEps = 10*test_precision<RealScalar>();
+  RealScalar eival_eps = (std::min)(test_precision<RealScalar>(),  NumTraits<Scalar>::dummy_precision()*20000);
   
   SelfAdjointEigenSolver<MatrixType> eiSymm(m);
   VERIFY_IS_EQUAL(eiSymm.info(), Success);
-  VERIFY((m.template selfadjointView<Lower>() * eiSymm.eigenvectors()).isApprox(
-          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
+  VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiSymm.eigenvectors(),
+                   eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal());
   VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiSymm.eigenvalues());
   VERIFY_IS_UNITARY(eiSymm.eigenvectors());
 
-  SelfAdjointEigenSolver<MatrixType> eiDirect;
-  eiDirect.computeDirect(m);  
-  VERIFY_IS_EQUAL(eiDirect.info(), Success);
-  VERIFY((m.template selfadjointView<Lower>() * eiDirect.eigenvectors()).isApprox(
-          eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal(), largerEps));
-  VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues());
-  VERIFY_IS_UNITARY(eiDirect.eigenvectors());
+  if(m.cols()<=4)
+  {
+    SelfAdjointEigenSolver<MatrixType> eiDirect;
+    eiDirect.computeDirect(m);  
+    VERIFY_IS_EQUAL(eiDirect.info(), Success);
+    VERIFY_IS_APPROX(eiSymm.eigenvalues(), eiDirect.eigenvalues());
+    if(! eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps) )
+    {
+      std::cerr << "reference eigenvalues: " << eiSymm.eigenvalues().transpose() << "\n"
+                << "obtained eigenvalues:  " << eiDirect.eigenvalues().transpose() << "\n"
+                << "diff:                  " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).transpose() << "\n"
+                << "error (eps):           " << (eiSymm.eigenvalues()-eiDirect.eigenvalues()).norm() / eiSymm.eigenvalues().norm() << "  (" << eival_eps << ")\n";
+    }
+    VERIFY(eiSymm.eigenvalues().isApprox(eiDirect.eigenvalues(), eival_eps));
+    VERIFY_IS_APPROX(m.template selfadjointView<Lower>() * eiDirect.eigenvectors(),
+                    eiDirect.eigenvectors() * eiDirect.eigenvalues().asDiagonal());
+    VERIFY_IS_APPROX(m.template selfadjointView<Lower>().eigenvalues(), eiDirect.eigenvalues());
+    VERIFY_IS_UNITARY(eiDirect.eigenvectors());
+  }
 }
 
 template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)