Move computation of eigenvalues from RealSchur to EigenSolver.

2 files changed, 33 insertions, 34 deletions

diff --git a/Eigen/src/Eigenvalues/EigenSolver.h b/Eigen/src/Eigenvalues/EigenSolver.h
index 44a0fd485..f8fc08efc 100644
--- a/Eigen/src/Eigenvalues/EigenSolver.h
+++ b/Eigen/src/Eigenvalues/EigenSolver.h
@@ -68,7 +68,9 @@
   * The documentation for EigenSolver(const MatrixType&) contains an example of
   * the typical use of this class.
   *
-  * \note this code was adapted from JAMA (public domain)
+  * \note The implementation is adapted from
+  * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
+  * Their code is based on EISPACK.
   *
   * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
   */
@@ -232,12 +234,13 @@ template<typename _MatrixType> class EigenSolver
       * The eigenvalues() and eigenvectors() functions can be used to retrieve
       * the computed eigendecomposition.
       *
-      * The matrix is first reduced to Schur form. The Schur decomposition is
-      * then used to compute the eigenvalues and eigenvectors.
+      * The matrix is first reduced to real Schur form using the RealSchur
+      * class. The Schur decomposition is then used to compute the eigenvalues
+      * and eigenvectors.
       *
       * The cost of the computation is dominated by the cost of the Schur
-      * decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ is the size of
-      * the matrix.
+      * decomposition, which is very approximately \f$ 25n^3 \f$ where 
+      * \f$ n \f$ is the size of the matrix.
       *
       * This method reuses of the allocated data in the EigenSolver object.
       *
@@ -311,12 +314,31 @@ EigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matr
 
   // Reduce to real Schur form.
   RealSchur<MatrixType> rs(matrix);
-  MatrixType matH = rs.matrixT();
+  MatrixType matT = rs.matrixT();
   m_eivec = rs.matrixU();
-  m_eivalues = rs.eigenvalues();
 
+  // Compute eigenvalues from matT
+  m_eivalues.resize(matrix.cols());
+  int i = 0;
+  while (i < matrix.cols()) 
+  {
+    if (i == matrix.cols() - 1 || matT.coeff(i+1, i) == Scalar(0)) 
+    {
+      m_eivalues.coeffRef(i) = matT.coeff(i, i);
+      ++i;
+    }
+    else
+    {
+      Scalar p = Scalar(0.5) * (matT.coeff(i, i) - matT.coeff(i+1, i+1));
+      Scalar z = ei_sqrt(ei_abs(p * p + matT.coeff(i+1, i) * matT.coeff(i, i+1)));
+      m_eivalues.coeffRef(i)   = ComplexScalar(matT.coeff(i+1, i+1) + p, z);
+      m_eivalues.coeffRef(i+1) = ComplexScalar(matT.coeff(i+1, i+1) + p, -z);
+      i += 2;
+    }
+  }
+  
   // Compute eigenvectors.
-  hqr2_step2(matH);
+  hqr2_step2(matT);
 
   m_isInitialized = true;
   return *this;
diff --git a/Eigen/src/Eigenvalues/RealSchur.h b/Eigen/src/Eigenvalues/RealSchur.h
index 17a3801ac..29569e802 100644
--- a/Eigen/src/Eigenvalues/RealSchur.h
+++ b/Eigen/src/Eigenvalues/RealSchur.h
@@ -53,7 +53,7 @@
   * given matrix. Alternatively, you can use the RealSchur(const MatrixType&)
   * constructor which computes the real Schur decomposition at construction
   * time. Once the decomposition is computed, you can use the matrixU() and
-  * matrixT() functions to retrieve the matrices U and V in the decomposition.
+  * matrixT() functions to retrieve the matrices U and T in the decomposition.
   *
   * The documentation of RealSchur(const MatrixType&) contains an example of
   * the typical use of this class.
@@ -93,7 +93,6 @@ template<typename _MatrixType> class RealSchur
     RealSchur(int size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime)
             : m_matT(size, size),
               m_matU(size, size), 
-              m_eivalues(size),
               m_isInitialized(false)
     { }
 
@@ -109,7 +108,6 @@ template<typename _MatrixType> class RealSchur
     RealSchur(const MatrixType& matrix)
             : m_matT(matrix.rows(),matrix.cols()),
               m_matU(matrix.rows(),matrix.cols()),
-              m_eivalues(matrix.rows()),
               m_isInitialized(false)
     {
       compute(matrix);
@@ -147,15 +145,6 @@ template<typename _MatrixType> class RealSchur
       return m_matT;
     }
   
-    /** \brief Returns vector of eigenvalues. 
-      *
-      * This function will likely be removed. */
-    const EigenvalueType& eigenvalues() const
-    {
-      ei_assert(m_isInitialized && "RealSchur is not initialized.");
-      return m_eivalues;
-    }
-  
     /** \brief Computes Schur decomposition of given matrix. 
       * 
       * \param[in]  matrix  Square matrix whose Schur decomposition is to be computed.
@@ -176,7 +165,6 @@ template<typename _MatrixType> class RealSchur
     
     MatrixType m_matT;
     MatrixType m_matU;
-    EigenvalueType m_eivalues;
     bool m_isInitialized;
 
     typedef Matrix<Scalar,3,1> Vector3s;
@@ -200,7 +188,6 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
   HessenbergDecomposition<MatrixType> hess(matrix);
   m_matT = hess.matrixH();
   m_matU = hess.matrixQ();
-  m_eivalues.resize(matrix.rows());
 
   // Step 2. Reduce to real Schur form  
   typedef Matrix<Scalar, ColsAtCompileTime, 1, Options, MaxColsAtCompileTime, 1> ColumnVectorType;
@@ -226,7 +213,6 @@ void RealSchur<MatrixType>::compute(const MatrixType& matrix)
       m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
       if (iu > 0) 
 	m_matT.coeffRef(iu, iu-1) = Scalar(0);
-      m_eivalues.coeffRef(iu) = ComplexScalar(m_matT.coeff(iu,iu), 0.0);
       iu--;
       iter = 0;
     }
@@ -289,15 +275,14 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
 
   // The eigenvalues of the 2x2 matrix [a b; c d] are 
   // trace +/- sqrt(discr/4) where discr = tr^2 - 4*det, tr = a + d, det = ad - bc
-  Scalar w = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
   Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
-  Scalar q = p * p + w;   // q = tr^2 / 4 - det = discr/4
-  Scalar z = ei_sqrt(ei_abs(q));
+  Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);   // q = tr^2 / 4 - det = discr/4
   m_matT.coeffRef(iu,iu) += exshift;
   m_matT.coeffRef(iu-1,iu-1) += exshift;
 
   if (q >= 0) // Two real eigenvalues
   {
+    Scalar z = ei_sqrt(ei_abs(q));
     PlanarRotation<Scalar> rot;
     if (p >= 0)
       rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
@@ -308,14 +293,6 @@ inline void RealSchur<MatrixType>::splitOffTwoRows(int iu, Scalar exshift)
     m_matT.block(0, 0, iu+1, size).applyOnTheRight(iu-1, iu, rot);
     m_matT.coeffRef(iu, iu-1) = Scalar(0); 
     m_matU.applyOnTheRight(iu-1, iu, rot);
-
-    m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu-1, iu-1), 0.0);
-    m_eivalues.coeffRef(iu)   = ComplexScalar(m_matT.coeff(iu, iu), 0.0);
-  }
-  else // // Pair of complex conjugate eigenvalues
-  {
-    m_eivalues.coeffRef(iu-1) = ComplexScalar(m_matT.coeff(iu,iu) + p, z);
-    m_eivalues.coeffRef(iu)   = ComplexScalar(m_matT.coeff(iu,iu) + p, -z);
   }
 
   if (iu > 1)