From a0cff1a2959208f8b767b8ff03b716e58ddb7ee1 Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Tue, 3 Jun 2008 18:03:55 +0000
Subject: [PATCH] fix eigenvectors computations :)

---
 Eigen/src/Core/Product.h              | 2 +-
 Eigen/src/QR/SelfAdjointEigenSolver.h | 7 +++----
 Eigen/src/QR/Tridiagonalization.h     | 4 ++++
 3 files changed, 8 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/Core/Product.h b/Eigen/src/Core/Product.h
index d29506dd3..3709ebf64 100644
--- a/Eigen/src/Core/Product.h
+++ b/Eigen/src/Core/Product.h
@@ -281,7 +281,7 @@ template<typename Lhs, typename Rhs, int EvalMode> class Product : ei_no_assignm
     {
       if ((Rhs::Flags&Diagonal)==Diagonal)
       {
-        assert(_LhsNested::Flags&RowMajorBit==0);
+        assert((_LhsNested::Flags&RowMajorBit)==0);
         return ei_pmul(m_lhs.template packetCoeff<LoadMode>(row, col), ei_pset1(m_rhs.coeff(col, col)));
       }
       else
diff --git a/Eigen/src/QR/SelfAdjointEigenSolver.h b/Eigen/src/QR/SelfAdjointEigenSolver.h
index 01b31e704..9b27a0d0c 100644
--- a/Eigen/src/QR/SelfAdjointEigenSolver.h
+++ b/Eigen/src/QR/SelfAdjointEigenSolver.h
@@ -45,7 +45,6 @@ template<typename _MatrixType> class SelfAdjointEigenSolver
     typedef typename MatrixType::Scalar Scalar;
     typedef typename NumTraits<Scalar>::Real RealScalar;
     typedef std::complex<RealScalar> Complex;
-//     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> EigenvalueType;
     typedef Matrix<RealScalar, MatrixType::ColsAtCompileTime, 1> RealVectorType;
     typedef Matrix<RealScalar, Dynamic, 1> RealVectorTypeX;
 
@@ -120,6 +119,9 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   m_eivalues.resize(n,1);
   m_eivec = matrix;
 
+  // FIXME, should tridiag be a local variable of this function or an attribute of SelfAdjointEigenSolver ?
+  // the latter avoids multiple memory allocation when the same SelfAdjointEigenSolver is used multiple times...
+  // (same for diag and subdiag)
   Tridiagonalization<MatrixType> tridiag(m_eivec);
   RealVectorType& diag = m_eivalues;
   RealVectorTypeX subdiag(n-1);
@@ -128,9 +130,6 @@ void SelfAdjointEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool
   if (computeEigenvectors)
     m_eivec = tridiag.matrixQ();
 
-  RealVectorTypeX gc(n);
-  RealVectorTypeX gs(n);
-
   int end = n-1;
   int start = 0;
   while (end>0)
diff --git a/Eigen/src/QR/Tridiagonalization.h b/Eigen/src/QR/Tridiagonalization.h
index 96a4292b3..6621d257b 100755
--- a/Eigen/src/QR/Tridiagonalization.h
+++ b/Eigen/src/QR/Tridiagonalization.h
@@ -197,6 +197,10 @@ void Tridiagonalization<MatrixType>::_compute(MatrixType& matA, CoeffVectorType&
     matA.col(i).coeffRef(i+1) = beta;
     hCoeffs.coeffRef(i) = (beta - v0) / beta;
   }
+  else
+  {
+    hCoeffs.coeffRef(n-2) = 0;
+  }
 }
 
 /** reconstructs and returns the matrix Q */