Main MRPT website > C++ reference
MRPT logo

EigenSolver.h

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



Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011