From 574ad9efbd80f91bd0db0cad94d42cd857d7dae7 Mon Sep 17 00:00:00 2001 From: Jitse Niesen Date: Mon, 12 Apr 2010 18:54:15 +0100 Subject: [PATCH] Move computation of eigenvalues from RealSchur to EigenSolver. --- Eigen/src/Eigenvalues/EigenSolver.h | 38 +++++++++++++++++++++++------ Eigen/src/Eigenvalues/RealSchur.h | 29 +++------------------- 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 + * JAMA (public domain). + * Their code is based on EISPACK. * * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver */ @@ -232,12 +234,13 @@ template 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& EigenSolver::compute(const MatrixType& matr // Reduce to real Schur form. RealSchur 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 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 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 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 class RealSchur MatrixType m_matT; MatrixType m_matU; - EigenvalueType m_eivalues; bool m_isInitialized; typedef Matrix Vector3s; @@ -200,7 +188,6 @@ void RealSchur::compute(const MatrixType& matrix) HessenbergDecomposition hess(matrix); m_matT = hess.matrixH(); m_matU = hess.matrixQ(); - m_eivalues.resize(matrix.rows()); // Step 2. Reduce to real Schur form typedef Matrix ColumnVectorType; @@ -226,7 +213,6 @@ void RealSchur::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::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 rot; if (p >= 0) rot.makeGivens(p + z, m_matT.coeff(iu, iu-1)); @@ -308,14 +293,6 @@ inline void RealSchur::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)