EigenSolver.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_EIGENSOLVER_H
00012 #define EIGEN_EIGENSOLVER_H
00013 
00014 #include "./RealSchur.h"
00015 
00016 namespace Eigen { 
00017 
00018 /** \eigenvalues_module \ingroup Eigenvalues_Module
00019   *
00020   *
00021   * \class EigenSolver
00022   *
00023   * \brief Computes eigenvalues and eigenvectors of general matrices
00024   *
00025   * \tparam _MatrixType the type of the matrix of which we are computing the
00026   * eigendecomposition; this is expected to be an instantiation of the Matrix
00027   * class template. Currently, only real matrices are supported.
00028   *
00029   * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars
00030   * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$.  If
00031   * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and
00032   * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V =
00033   * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we
00034   * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition.
00035   *
00036   * The eigenvalues and eigenvectors of a matrix may be complex, even when the
00037   * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D
00038   * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the
00039   * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to
00040   * have blocks of the form
00041   * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f]
00042   * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal.  These
00043   * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call
00044   * this variant of the eigendecomposition the pseudo-eigendecomposition.
00045   *
00046   * Call the function compute() to compute the eigenvalues and eigenvectors of
00047   * a given matrix. Alternatively, you can use the 
00048   * EigenSolver(const MatrixType&, bool) constructor which computes the
00049   * eigenvalues and eigenvectors at construction time. Once the eigenvalue and
00050   * eigenvectors are computed, they can be retrieved with the eigenvalues() and
00051   * eigenvectors() functions. The pseudoEigenvalueMatrix() and
00052   * pseudoEigenvectors() methods allow the construction of the
00053   * pseudo-eigendecomposition.
00054   *
00055   * The documentation for EigenSolver(const MatrixType&, bool) contains an
00056   * example of the typical use of this class.
00057   *
00058   * \note The implementation is adapted from
00059   * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain).
00060   * Their code is based on EISPACK.
00061   *
00062   * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver
00063   */
00064 template<typename _MatrixType> class EigenSolver 
00065 {
00066   public:
00067 
00068     /** \brief Synonym for the template parameter \p _MatrixType. */
00069     typedef _MatrixType MatrixType;
00070 
00071     enum {
00072       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00073       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00074       Options = MatrixType::Options,
00075       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00076       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00077     };
00078 
00079     /** \brief Scalar type for matrices of type #MatrixType. */
00080     typedef typename MatrixType::Scalar Scalar;
00081     typedef typename NumTraits<Scalar>::Real RealScalar;
00082     typedef typename MatrixType::Index Index;
00083 
00084     /** \brief Complex scalar type for #MatrixType. 
00085       *
00086       * This is \c std::complex<Scalar> if #Scalar is real (e.g.,
00087       * \c float or \c double) and just \c Scalar if #Scalar is
00088       * complex.
00089       */
00090     typedef std::complex<RealScalar> ComplexScalar;
00091 
00092     /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 
00093       *
00094       * This is a column vector with entries of type #ComplexScalar.
00095       * The length of the vector is the size of #MatrixType.
00096       */
00097     typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1>  EigenvalueType;
00098 
00099     /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 
00100       *
00101       * This is a square matrix with entries of type #ComplexScalar. 
00102       * The size is the same as the size of #MatrixType.
00103       */
00104     typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime>  EigenvectorsType;
00105 
00106     /** \brief Default constructor.
00107       *
00108       * The default constructor is useful in cases in which the user intends to
00109       * perform decompositions via EigenSolver::compute(const MatrixType&, bool).
00110       *
00111       * \sa compute() for an example.
00112       */
00113  EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {}
00114 
00115     /** \brief Default constructor with memory preallocation
00116       *
00117       * Like the default constructor but with preallocation of the internal data
00118       * according to the specified problem \a size.
00119       * \sa EigenSolver()
00120       */
00121     EigenSolver(Index size)
00122       : m_eivec(size, size),
00123         m_eivalues(size),
00124         m_isInitialized(false),
00125         m_eigenvectorsOk(false),
00126         m_realSchur(size),
00127         m_matT(size, size), 
00128         m_tmp(size)
00129     {}
00130 
00131     /** \brief Constructor; computes eigendecomposition of given matrix. 
00132       * 
00133       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00134       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00135       *    eigenvalues are computed; if false, only the eigenvalues are
00136       *    computed. 
00137       *
00138       * This constructor calls compute() to compute the eigenvalues
00139       * and eigenvectors.
00140       *
00141       * Example: \include EigenSolver_EigenSolver_MatrixType.cpp
00142       * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out
00143       *
00144       * \sa compute()
00145       */
00146     EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true)
00147       : m_eivec(matrix.rows(), matrix.cols()),
00148         m_eivalues(matrix.cols()),
00149         m_isInitialized(false),
00150         m_eigenvectorsOk(false),
00151         m_realSchur(matrix.cols()),
00152         m_matT(matrix.rows(), matrix.cols()), 
00153         m_tmp(matrix.cols())
00154     {
00155       compute(matrix, computeEigenvectors);
00156     }
00157 
00158     /** \brief Returns the eigenvectors of given matrix. 
00159       *
00160       * \returns  %Matrix whose columns are the (possibly complex) eigenvectors.
00161       *
00162       * \pre Either the constructor 
00163       * EigenSolver(const MatrixType&,bool) or the member function
00164       * compute(const MatrixType&, bool) has been called before, and
00165       * \p computeEigenvectors was set to true (the default).
00166       *
00167       * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
00168       * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
00169       * eigenvectors are normalized to have (Euclidean) norm equal to one. The
00170       * matrix returned by this function is the matrix \f$ V \f$ in the
00171       * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists.
00172       *
00173       * Example: \include EigenSolver_eigenvectors.cpp
00174       * Output: \verbinclude EigenSolver_eigenvectors.out
00175       *
00176       * \sa eigenvalues(), pseudoEigenvectors()
00177       */
00178     EigenvectorsType eigenvectors() const;
00179 
00180     /** \brief Returns the pseudo-eigenvectors of given matrix. 
00181       *
00182       * \returns  Const reference to matrix whose columns are the pseudo-eigenvectors.
00183       *
00184       * \pre Either the constructor 
00185       * EigenSolver(const MatrixType&,bool) or the member function
00186       * compute(const MatrixType&, bool) has been called before, and
00187       * \p computeEigenvectors was set to true (the default).
00188       *
00189       * The real matrix \f$ V \f$ returned by this function and the
00190       * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix()
00191       * satisfy \f$ AV = VD \f$.
00192       *
00193       * Example: \include EigenSolver_pseudoEigenvectors.cpp
00194       * Output: \verbinclude EigenSolver_pseudoEigenvectors.out
00195       *
00196       * \sa pseudoEigenvalueMatrix(), eigenvectors()
00197       */
00198     const MatrixType& pseudoEigenvectors() const
00199     {
00200       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00201       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00202       return m_eivec;
00203     }
00204 
00205     /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition.
00206       *
00207       * \returns  A block-diagonal matrix.
00208       *
00209       * \pre Either the constructor 
00210       * EigenSolver(const MatrixType&,bool) or the member function
00211       * compute(const MatrixType&, bool) has been called before.
00212       *
00213       * The matrix \f$ D \f$ returned by this function is real and
00214       * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2
00215       * blocks of the form
00216       * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$.
00217       * These blocks are not sorted in any particular order.
00218       * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by
00219       * pseudoEigenvectors() satisfy \f$ AV = VD \f$.
00220       *
00221       * \sa pseudoEigenvectors() for an example, eigenvalues()
00222       */
00223     MatrixType pseudoEigenvalueMatrix() const;
00224 
00225     /** \brief Returns the eigenvalues of given matrix. 
00226       *
00227       * \returns A const reference to the column vector containing the eigenvalues.
00228       *
00229       * \pre Either the constructor 
00230       * EigenSolver(const MatrixType&,bool) or the member function
00231       * compute(const MatrixType&, bool) has been called before.
00232       *
00233       * The eigenvalues are repeated according to their algebraic multiplicity,
00234       * so there are as many eigenvalues as rows in the matrix. The eigenvalues 
00235       * are not sorted in any particular order.
00236       *
00237       * Example: \include EigenSolver_eigenvalues.cpp
00238       * Output: \verbinclude EigenSolver_eigenvalues.out
00239       *
00240       * \sa eigenvectors(), pseudoEigenvalueMatrix(),
00241       *     MatrixBase::eigenvalues()
00242       */
00243     const EigenvalueType & eigenvalues() const
00244     {
00245       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00246       return m_eivalues;
00247     }
00248 
00249     /** \brief Computes eigendecomposition of given matrix. 
00250       * 
00251       * \param[in]  matrix  Square matrix whose eigendecomposition is to be computed.
00252       * \param[in]  computeEigenvectors  If true, both the eigenvectors and the
00253       *    eigenvalues are computed; if false, only the eigenvalues are
00254       *    computed. 
00255       * \returns    Reference to \c *this
00256       *
00257       * This function computes the eigenvalues of the real matrix \p matrix.
00258       * The eigenvalues() function can be used to retrieve them.  If 
00259       * \p computeEigenvectors is true, then the eigenvectors are also computed
00260       * and can be retrieved by calling eigenvectors().
00261       *
00262       * The matrix is first reduced to real Schur form using the RealSchur
00263       * class. The Schur decomposition is then used to compute the eigenvalues
00264       * and eigenvectors.
00265       *
00266       * The cost of the computation is dominated by the cost of the
00267       * Schur decomposition, which is very approximately \f$ 25n^3 \f$
00268       * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors 
00269       * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false.
00270       *
00271       * This method reuses of the allocated data in the EigenSolver object.
00272       *
00273       * Example: \include EigenSolver_compute.cpp
00274       * Output: \verbinclude EigenSolver_compute.out
00275       */
00276     EigenSolver & compute(const MatrixType& matrix, bool computeEigenvectors = true);
00277 
00278     ComputationInfo info() const
00279     {
00280       eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00281       return m_realSchur.info();
00282     }
00283 
00284     /** \brief Sets the maximum number of iterations allowed. */
00285     EigenSolver & setMaxIterations(Index maxIters)
00286     {
00287       m_realSchur.setMaxIterations(maxIters);
00288       return *this;
00289     }
00290 
00291     /** \brief Returns the maximum number of iterations. */
00292     Index getMaxIterations()
00293     {
00294       return m_realSchur.getMaxIterations();
00295     }
00296 
00297   private:
00298     void doComputeEigenvectors();
00299 
00300   protected:
00301     
00302     static void check_template_parameters()
00303     {
00304       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00305       EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL);
00306     }
00307     
00308     MatrixType m_eivec;
00309     EigenvalueType m_eivalues;
00310     bool m_isInitialized;
00311     bool m_eigenvectorsOk;
00312     RealSchur<MatrixType> m_realSchur;
00313     MatrixType m_matT;
00314 
00315     typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
00316     ColumnVectorType m_tmp;
00317 };
00318 
00319 template<typename MatrixType>
00320 MatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const
00321 {
00322   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00323   Index n = m_eivalues.rows();
00324   MatrixType matD = MatrixType::Zero(n,n);
00325   for (Index i=0; i<n; ++i)
00326   {
00327     if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i))))
00328       matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i));
00329     else
00330     {
00331       matD.template block<2,2>(i,i) <<  numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)),
00332                                        -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i));
00333       ++i;
00334     }
00335   }
00336   return matD;
00337 }
00338 
00339 template<typename MatrixType>
00340 typename EigenSolver<MatrixType>::EigenvectorsType  EigenSolver<MatrixType>::eigenvectors() const
00341 {
00342   eigen_assert(m_isInitialized && "EigenSolver is not initialized.");
00343   eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
00344   Index n = m_eivec.cols();
00345   EigenvectorsType  matV(n,n);
00346   for (Index j=0; j<n; ++j)
00347   {
00348     if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n)
00349     {
00350       // we have a real eigen value
00351       matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>();
00352       matV.col(j).normalize();
00353     }
00354     else
00355     {
00356       // we have a pair of complex eigen values
00357       for (Index i=0; i<n; ++i)
00358       {
00359         matV.coeffRef(i,j)   = ComplexScalar(m_eivec.coeff(i,j),  m_eivec.coeff(i,j+1));
00360         matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1));
00361       }
00362       matV.col(j).normalize();
00363       matV.col(j+1).normalize();
00364       ++j;
00365     }
00366   }
00367   return matV;
00368 }
00369 
00370 template<typename MatrixType>
00371 EigenSolver<MatrixType> & 
00372 EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors)
00373 {
00374   check_template_parameters();
00375   
00376   using std::sqrt;
00377   using std::abs;
00378   eigen_assert(matrix.cols() == matrix.rows());
00379 
00380   // Reduce to real Schur form.
00381   m_realSchur.compute(matrix, computeEigenvectors);
00382 
00383   if (m_realSchur.info() == Success)
00384   {
00385     m_matT = m_realSchur.matrixT();
00386     if (computeEigenvectors)
00387       m_eivec = m_realSchur.matrixU();
00388   
00389     // Compute eigenvalues from matT
00390     m_eivalues.resize(matrix.cols());
00391     Index i = 0;
00392     while (i < matrix.cols()) 
00393     {
00394       if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 
00395       {
00396         m_eivalues.coeffRef(i) = m_matT.coeff(i, i);
00397         ++i;
00398       }
00399       else
00400       {
00401         Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1));
00402         Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1)));
00403         m_eivalues.coeffRef(i)   = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z);
00404         m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z);
00405         i += 2;
00406       }
00407     }
00408     
00409     // Compute eigenvectors.
00410     if (computeEigenvectors)
00411       doComputeEigenvectors();
00412   }
00413 
00414   m_isInitialized = true;
00415   m_eigenvectorsOk = computeEigenvectors;
00416 
00417   return *this;
00418 }
00419 
00420 // Complex scalar division.
00421 template<typename Scalar>
00422 std::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi)
00423 {
00424   using std::abs;
00425   Scalar r,d;
00426   if (abs(yr) > abs(yi))
00427   {
00428       r = yi/yr;
00429       d = yr + r*yi;
00430       return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d);
00431   }
00432   else
00433   {
00434       r = yr/yi;
00435       d = yi + r*yr;
00436       return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d);
00437   }
00438 }
00439 
00440 
00441 template<typename MatrixType>
00442 void EigenSolver<MatrixType>::doComputeEigenvectors()
00443 {
00444   using std::abs;
00445   const Index size = m_eivec.cols();
00446   const Scalar eps = NumTraits<Scalar>::epsilon();
00447 
00448   // inefficient! this is already computed in RealSchur
00449   Scalar norm(0);
00450   for (Index j = 0; j < size; ++j)
00451   {
00452     norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
00453   }
00454   
00455   // Backsubstitute to find vectors of upper triangular form
00456   if (norm == 0.0)
00457   {
00458     return;
00459   }
00460 
00461   for (Index n = size-1; n >= 0; n--)
00462   {
00463     Scalar p = m_eivalues.coeff(n).real();
00464     Scalar q = m_eivalues.coeff(n).imag();
00465 
00466     // Scalar vector
00467     if (q == Scalar(0))
00468     {
00469       Scalar lastr(0), lastw(0);
00470       Index l = n;
00471 
00472       m_matT.coeffRef(n,n) = 1.0;
00473       for (Index i = n-1; i >= 0; i--)
00474       {
00475         Scalar w = m_matT.coeff(i,i) - p;
00476         Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00477 
00478         if (m_eivalues.coeff(i).imag() < 0.0)
00479         {
00480           lastw = w;
00481           lastr = r;
00482         }
00483         else
00484         {
00485           l = i;
00486           if (m_eivalues.coeff(i).imag() == 0.0)
00487           {
00488             if (w != 0.0)
00489               m_matT.coeffRef(i,n) = -r / w;
00490             else
00491               m_matT.coeffRef(i,n) = -r / (eps * norm);
00492           }
00493           else // Solve real equations
00494           {
00495             Scalar x = m_matT.coeff(i,i+1);
00496             Scalar y = m_matT.coeff(i+1,i);
00497             Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag();
00498             Scalar t = (x * lastr - lastw * r) / denom;
00499             m_matT.coeffRef(i,n) = t;
00500             if (abs(x) > abs(lastw))
00501               m_matT.coeffRef(i+1,n) = (-r - w * t) / x;
00502             else
00503               m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw;
00504           }
00505 
00506           // Overflow control
00507           Scalar t = abs(m_matT.coeff(i,n));
00508           if ((eps * t) * t > Scalar(1))
00509             m_matT.col(n).tail(size-i) /= t;
00510         }
00511       }
00512     }
00513     else if (q < Scalar(0) && n > 0) // Complex vector
00514     {
00515       Scalar lastra(0), lastsa(0), lastw(0);
00516       Index l = n-1;
00517 
00518       // Last vector component imaginary so matrix is triangular
00519       if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n)))
00520       {
00521         m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1);
00522         m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1);
00523       }
00524       else
00525       {
00526         std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q);
00527         m_matT.coeffRef(n-1,n-1) = numext::real(cc);
00528         m_matT.coeffRef(n-1,n) = numext::imag(cc);
00529       }
00530       m_matT.coeffRef(n,n-1) = 0.0;
00531       m_matT.coeffRef(n,n) = 1.0;
00532       for (Index i = n-2; i >= 0; i--)
00533       {
00534         Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1));
00535         Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1));
00536         Scalar w = m_matT.coeff(i,i) - p;
00537 
00538         if (m_eivalues.coeff(i).imag() < 0.0)
00539         {
00540           lastw = w;
00541           lastra = ra;
00542           lastsa = sa;
00543         }
00544         else
00545         {
00546           l = i;
00547           if (m_eivalues.coeff(i).imag() == RealScalar(0))
00548           {
00549             std::complex<Scalar> cc = cdiv(-ra,-sa,w,q);
00550             m_matT.coeffRef(i,n-1) = numext::real(cc);
00551             m_matT.coeffRef(i,n) = numext::imag(cc);
00552           }
00553           else
00554           {
00555             // Solve complex equations
00556             Scalar x = m_matT.coeff(i,i+1);
00557             Scalar y = m_matT.coeff(i+1,i);
00558             Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q;
00559             Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q;
00560             if ((vr == 0.0) && (vi == 0.0))
00561               vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw));
00562 
00563             std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi);
00564             m_matT.coeffRef(i,n-1) = numext::real(cc);
00565             m_matT.coeffRef(i,n) = numext::imag(cc);
00566             if (abs(x) > (abs(lastw) + abs(q)))
00567             {
00568               m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x;
00569               m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x;
00570             }
00571             else
00572             {
00573               cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q);
00574               m_matT.coeffRef(i+1,n-1) = numext::real(cc);
00575               m_matT.coeffRef(i+1,n) = numext::imag(cc);
00576             }
00577           }
00578 
00579           // Overflow control
00580           using std::max;
00581           Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n)));
00582           if ((eps * t) * t > Scalar(1))
00583             m_matT.block(i, n-1, size-i, 2) /= t;
00584 
00585         }
00586       }
00587       
00588       // We handled a pair of complex conjugate eigenvalues, so need to skip them both
00589       n--;
00590     }
00591     else
00592     {
00593       eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen
00594     }
00595   }
00596 
00597   // Back transformation to get eigenvectors of original matrix
00598   for (Index j = size-1; j >= 0; j--)
00599   {
00600     m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1);
00601     m_eivec.col(j) = m_tmp;
00602   }
00603 }
00604 
00605 } // end namespace Eigen
00606 
00607 #endif // EIGEN_EIGENSOLVER_H