From 374deaed5fcebe08922e2451107a9727851b3e05 Mon Sep 17 00:00:00 2001
From: Benoit Jacob <jacob.benoit.1@gmail.com>
Date: Mon, 31 Jan 2011 08:36:14 -0500
Subject: [PATCH] make eigen2 eigensolver test pass

---
 .../src/Eigenvalues/SelfAdjointEigenSolver.h  | 33 +++++++++++++++++++
 test/eigen2/eigen2_eigensolver.cpp            | 10 +++---
 2 files changed, 38 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
index 83d9c015c..3f3b08193 100644
--- a/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
+++ b/Eigen/src/Eigenvalues/SelfAdjointEigenSolver.h
@@ -29,6 +29,9 @@
 #include "./EigenvaluesCommon.h"
 #include "./Tridiagonalization.h"
 
+template<typename _MatrixType>
+class GeneralizedSelfAdjointEigenSolver;
+
 /** \eigenvalues_module \ingroup Eigenvalues_Module
   *
   *
@@ -310,6 +313,36 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
       */
     static const int m_maxIterations = 30;
 
+    #ifdef EIGEN2_SUPPORT
+    SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors)
+      : m_eivec(matrix.rows(), matrix.cols()),
+        m_eivalues(matrix.cols()),
+        m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
+        m_isInitialized(false)
+    {
+      compute(matrix, computeEigenvectors);
+    }
+    
+    SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
+        : m_eivec(matA.cols(), matA.cols()),
+          m_eivalues(matA.cols()),
+          m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1),
+          m_isInitialized(false)
+    {
+      static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
+    }
+    
+    void compute(const MatrixType& matrix, bool computeEigenvectors)
+    {
+      compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
+    }
+
+    void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true)
+    {
+      compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly);
+    }
+    #endif // EIGEN2_SUPPORT
+
   protected:
     MatrixType m_eivec;
     RealVectorType m_eivalues;
diff --git a/test/eigen2/eigen2_eigensolver.cpp b/test/eigen2/eigen2_eigensolver.cpp
index 17111d559..e9b64cc4b 100644
--- a/test/eigen2/eigen2_eigensolver.cpp
+++ b/test/eigen2/eigen2_eigensolver.cpp
@@ -75,7 +75,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     convert(gEvec, _evec);
 
     // test gsl itself !
-    VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal().eval(), largerEps));
+    VERIFY((symmA * _evec).isApprox(_evec * _eval.asDiagonal(), largerEps));
 
     // compare with eigen
     VERIFY_IS_APPROX(_eval, eiSymm.eigenvalues());
@@ -86,7 +86,7 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
     convert(gEval, _eval);
     convert(gEvec, _evec);
     // test GSL itself:
-    VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal().eval()), largerEps));
+    VERIFY((symmA * _evec).isApprox(symmB * (_evec * _eval.asDiagonal()), largerEps));
 
     // compare with eigen
 //     std::cerr << _eval.transpose() << "\n" << eiSymmGen.eigenvalues().transpose() << "\n\n";
@@ -102,11 +102,11 @@ template<typename MatrixType> void selfadjointeigensolver(const MatrixType& m)
   #endif
 
   VERIFY((symmA * eiSymm.eigenvectors()).isApprox(
-          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal().eval(), largerEps));
+          eiSymm.eigenvectors() * eiSymm.eigenvalues().asDiagonal(), largerEps));
 
   // generalized eigen problem Ax = lBx
   VERIFY((symmA * eiSymmGen.eigenvectors()).isApprox(
-          symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal().eval()), largerEps));
+          symmB * (eiSymmGen.eigenvectors() * eiSymmGen.eigenvalues().asDiagonal()), largerEps));
 
   MatrixType sqrtSymmA = eiSymm.operatorSqrt();
   VERIFY_IS_APPROX(symmA, sqrtSymmA*sqrtSymmA);
@@ -141,7 +141,7 @@ template<typename MatrixType> void eigensolver(const MatrixType& m)
   EigenSolver<MatrixType> ei1(a);
   VERIFY_IS_APPROX(a * ei1.pseudoEigenvectors(), ei1.pseudoEigenvectors() * ei1.pseudoEigenvalueMatrix());
   VERIFY_IS_APPROX(a.template cast<Complex>() * ei1.eigenvectors(),
-                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal().eval());
+                   ei1.eigenvectors() * ei1.eigenvalues().asDiagonal());
 
 }