src/Eigenvalues/SelfAdjointEigenSolver.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H
0012 #define EIGEN_SELFADJOINTEIGENSOLVER_H
0013
0014 #include "./Tridiagonalization.h"
0015
0016 namespace Eigen {
0017
0018 template<typename _MatrixType>
0019 class GeneralizedSelfAdjointEigenSolver;
0020
0021 namespace internal {
0022 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues;
0023
0024 template<typename MatrixType, typename DiagType, typename SubDiagType>
0025 EIGEN_DEVICE_FUNC
0026 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec);
0027 }
0028
0029 /** \eigenvalues_module \ingroup Eigenvalues_Module
0030   *
0031   *
0032   * \class SelfAdjointEigenSolver
0033   *
0034   * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices
0035   *
0036   * \tparam _MatrixType the type of the matrix of which we are computing the
0037   * eigendecomposition; this is expected to be an instantiation of the Matrix
0038   * class template.
0039   *
0040   * A matrix \f$ A \f$ is selfadjoint if it equals its adjoint. For real
0041   * matrices, this means that the matrix is symmetric: it equals its
0042   * transpose. This class computes the eigenvalues and eigenvectors of a
0043   * selfadjoint matrix. These are the scalars \f$ \lambda \f$ and vectors
0044   * \f$ v \f$ such that \f$ Av = \lambda v \f$.  The eigenvalues of a
0045   * selfadjoint matrix are always real. If \f$ D \f$ is a diagonal matrix with
0046   * the eigenvalues on the diagonal, and \f$ V \f$ is a matrix with the
0047   * eigenvectors as its columns, then \f$ A = V D V^{-1} \f$. This is called the
0048   * eigendecomposition.
0049   *
0050   * For a selfadjoint matrix, \f$ V \f$ is unitary, meaning its inverse is equal
0051   * to its adjoint, \f$ V^{-1} = V^{\dagger} \f$. If \f$ A \f$ is real, then
0052   * \f$ V \f$ is also real and therefore orthogonal, meaning its inverse is
0053   * equal to its transpose, \f$ V^{-1} = V^T \f$.
0054   *
0055   * The algorithm exploits the fact that the matrix is selfadjoint, making it
0056   * faster and more accurate than the general purpose eigenvalue algorithms
0057   * implemented in EigenSolver and ComplexEigenSolver.
0058   *
0059   * Only the \b lower \b triangular \b part of the input matrix is referenced.
0060   *
0061   * Call the function compute() to compute the eigenvalues and eigenvectors of
0062   * a given matrix. Alternatively, you can use the
0063   * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes
0064   * the eigenvalues and eigenvectors at construction time. Once the eigenvalue
0065   * and eigenvectors are computed, they can be retrieved with the eigenvalues()
0066   * and eigenvectors() functions.
0067   *
0068   * The documentation for SelfAdjointEigenSolver(const MatrixType&, int)
0069   * contains an example of the typical use of this class.
0070   *
0071   * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and
0072   * the likes, see the class GeneralizedSelfAdjointEigenSolver.
0073   *
0074   * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver
0075   */
0076 template<typename _MatrixType> class SelfAdjointEigenSolver
0077 {
0078   public:
0079
0080     typedef _MatrixType MatrixType;
0081     enum {
0082       Size = MatrixType::RowsAtCompileTime,
0083       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0084       Options = MatrixType::Options,
0085       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0086     };
0087
0088     /** \brief Scalar type for matrices of type \p _MatrixType. */
0089     typedef typename MatrixType::Scalar Scalar;
0090     typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
0091
0092     typedef Matrix<Scalar,Size,Size,ColMajor,MaxColsAtCompileTime,MaxColsAtCompileTime> EigenvectorsType;
0093
0094     /** \brief Real scalar type for \p _MatrixType.
0095       *
0096       * This is just \c Scalar if #Scalar is real (e.g., \c float or
0097       * \c double), and the type of the real part of \c Scalar if #Scalar is
0098       * complex.
0099       */
0100     typedef typename NumTraits<Scalar>::Real RealScalar;
0101
0102     friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>;
0103
0104     /** \brief Type for vector of eigenvalues as returned by eigenvalues().
0105       *
0106       * This is a column vector with entries of type #RealScalar.
0107       * The length of the vector is the size of \p _MatrixType.
0108       */
0109     typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType;
0110     typedef Tridiagonalization<MatrixType> TridiagonalizationType;
0111     typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType;
0112
0113     /** \brief Default constructor for fixed-size matrices.
0114       *
0115       * The default constructor is useful in cases in which the user intends to
0116       * perform decompositions via compute(). This constructor
0117       * can only be used if \p _MatrixType is a fixed-size matrix; use
0118       * SelfAdjointEigenSolver(Index) for dynamic-size matrices.
0119       *
0120       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp
0121       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out
0122       */
0123     EIGEN_DEVICE_FUNC
0124     SelfAdjointEigenSolver()
0125         : m_eivec(),
0126           m_eivalues(),
0127           m_subdiag(),
0128           m_hcoeffs(),
0129           m_info(InvalidInput),
0130           m_isInitialized(false),
0131           m_eigenvectorsOk(false)
0132     { }
0133
0134     /** \brief Constructor, pre-allocates memory for dynamic-size matrices.
0135       *
0136       * \param [in]  size  Positive integer, size of the matrix whose
0137       * eigenvalues and eigenvectors will be computed.
0138       *
0139       * This constructor is useful for dynamic-size matrices, when the user
0140       * intends to perform decompositions via compute(). The \p size
0141       * parameter is only used as a hint. It is not an error to give a wrong
0142       * \p size, but it may impair performance.
0143       *
0144       * \sa compute() for an example
0145       */
0146     EIGEN_DEVICE_FUNC
0147     explicit SelfAdjointEigenSolver(Index size)
0148         : m_eivec(size, size),
0149           m_eivalues(size),
0150           m_subdiag(size > 1 ? size - 1 : 1),
0151           m_hcoeffs(size > 1 ? size - 1 : 1),
0152           m_isInitialized(false),
0153           m_eigenvectorsOk(false)
0154     {}
0155
0156     /** \brief Constructor; computes eigendecomposition of given matrix.
0157       *
0158       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
0159       *    be computed. Only the lower triangular part of the matrix is referenced.
0160       * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0161       *
0162       * This constructor calls compute(const MatrixType&, int) to compute the
0163       * eigenvalues of the matrix \p matrix. The eigenvectors are computed if
0164       * \p options equals #ComputeEigenvectors.
0165       *
0166       * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp
0167       * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out
0168       *
0169       * \sa compute(const MatrixType&, int)
0170       */
0171     template<typename InputType>
0172     EIGEN_DEVICE_FUNC
0173     explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors)
0174       : m_eivec(matrix.rows(), matrix.cols()),
0175         m_eivalues(matrix.cols()),
0176         m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1),
0177         m_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1),
0178         m_isInitialized(false),
0179         m_eigenvectorsOk(false)
0180     {
0181       compute(matrix.derived(), options);
0182     }
0183
0184     /** \brief Computes eigendecomposition of given matrix.
0185       *
0186       * \param[in]  matrix  Selfadjoint matrix whose eigendecomposition is to
0187       *    be computed. Only the lower triangular part of the matrix is referenced.
0188       * \param[in]  options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0189       * \returns    Reference to \c *this
0190       *
0191       * This function computes the eigenvalues of \p matrix.  The eigenvalues()
0192       * function can be used to retrieve them.  If \p options equals #ComputeEigenvectors,
0193       * then the eigenvectors are also computed and can be retrieved by
0194       * calling eigenvectors().
0195       *
0196       * This implementation uses a symmetric QR algorithm. The matrix is first
0197       * reduced to tridiagonal form using the Tridiagonalization class. The
0198       * tridiagonal matrix is then brought to diagonal form with implicit
0199       * symmetric QR steps with Wilkinson shift. Details can be found in
0200       * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>.
0201       *
0202       * The cost of the computation is about \f$ 9n^3 \f$ if the eigenvectors
0203       * are required and \f$ 4n^3/3 \f$ if they are not required.
0204       *
0205       * This method reuses the memory in the SelfAdjointEigenSolver object that
0206       * was allocated when the object was constructed, if the size of the
0207       * matrix does not change.
0208       *
0209       * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp
0210       * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out
0211       *
0212       * \sa SelfAdjointEigenSolver(const MatrixType&, int)
0213       */
0214     template<typename InputType>
0215     EIGEN_DEVICE_FUNC
0216     SelfAdjointEigenSolver& compute(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors);
0217
0218     /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm
0219       *
0220       * This is a variant of compute(const MatrixType&, int options) which
0221       * directly solves the underlying polynomial equation.
0222       *
0223       * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d).
0224       *
0225       * This method is usually significantly faster than the QR iterative algorithm
0226       * but it might also be less accurate. It is also worth noting that
0227       * for 3x3 matrices it involves trigonometric operations which are
0228       * not necessarily available for all scalar types.
0229       *
0230       * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues:
0231       *   - double: 1e-8
0232       *   - float:  1e-3
0233       *
0234       * \sa compute(const MatrixType&, int options)
0235       */
0236     EIGEN_DEVICE_FUNC
0237     SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors);
0238
0239     /**
0240       *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix
0241       *
0242       * \param[in] diag The vector containing the diagonal of the matrix.
0243       * \param[in] subdiag The subdiagonal of the matrix.
0244       * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly.
0245       * \returns Reference to \c *this
0246       *
0247       * This function assumes that the matrix has been reduced to tridiagonal form.
0248       *
0249       * \sa compute(const MatrixType&, int) for more information
0250       */
0251     SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors);
0252
0253     /** \brief Returns the eigenvectors of given matrix.
0254       *
0255       * \returns  A const reference to the matrix whose columns are the eigenvectors.
0256       *
0257       * \pre The eigenvectors have been computed before.
0258       *
0259       * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding
0260       * to eigenvalue number \f$ k \f$ as returned by eigenvalues().  The
0261       * eigenvectors are normalized to have (Euclidean) norm equal to one. If
0262       * this object was used to solve the eigenproblem for the selfadjoint
0263       * matrix \f$ A \f$, then the matrix returned by this function is the
0264       * matrix \f$ V \f$ in the eigendecomposition \f$ A = V D V^{-1} \f$.
0265       *
0266       * For a selfadjoint matrix, \f$ V \f$ is unitary, meaning its inverse is equal
0267       * to its adjoint, \f$ V^{-1} = V^{\dagger} \f$. If \f$ A \f$ is real, then
0268       * \f$ V \f$ is also real and therefore orthogonal, meaning its inverse is
0269       * equal to its transpose, \f$ V^{-1} = V^T \f$.
0270       *
0271       * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp
0272       * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out
0273       *
0274       * \sa eigenvalues()
0275       */
0276     EIGEN_DEVICE_FUNC
0277     const EigenvectorsType& eigenvectors() const
0278     {
0279       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0280       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0281       return m_eivec;
0282     }
0283
0284     /** \brief Returns the eigenvalues of given matrix.
0285       *
0286       * \returns A const reference to the column vector containing the eigenvalues.
0287       *
0288       * \pre The eigenvalues have been computed before.
0289       *
0290       * The eigenvalues are repeated according to their algebraic multiplicity,
0291       * so there are as many eigenvalues as rows in the matrix. The eigenvalues
0292       * are sorted in increasing order.
0293       *
0294       * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp
0295       * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out
0296       *
0297       * \sa eigenvectors(), MatrixBase::eigenvalues()
0298       */
0299     EIGEN_DEVICE_FUNC
0300     const RealVectorType& eigenvalues() const
0301     {
0302       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0303       return m_eivalues;
0304     }
0305
0306     /** \brief Computes the positive-definite square root of the matrix.
0307       *
0308       * \returns the positive-definite square root of the matrix
0309       *
0310       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
0311       * have been computed before.
0312       *
0313       * The square root of a positive-definite matrix \f$ A \f$ is the
0314       * positive-definite matrix whose square equals \f$ A \f$. This function
0315       * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the
0316       * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$.
0317       *
0318       * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp
0319       * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out
0320       *
0321       * \sa operatorInverseSqrt(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a>
0322       */
0323     EIGEN_DEVICE_FUNC
0324     MatrixType operatorSqrt() const
0325     {
0326       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0327       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0328       return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint();
0329     }
0330
0331     /** \brief Computes the inverse square root of the matrix.
0332       *
0333       * \returns the inverse positive-definite square root of the matrix
0334       *
0335       * \pre The eigenvalues and eigenvectors of a positive-definite matrix
0336       * have been computed before.
0337       *
0338       * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to
0339       * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is
0340       * cheaper than first computing the square root with operatorSqrt() and
0341       * then its inverse with MatrixBase::inverse().
0342       *
0343       * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp
0344       * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out
0345       *
0346       * \sa operatorSqrt(), MatrixBase::inverse(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a>
0347       */
0348     EIGEN_DEVICE_FUNC
0349     MatrixType operatorInverseSqrt() const
0350     {
0351       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0352       eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues.");
0353       return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint();
0354     }
0355
0356     /** \brief Reports whether previous computation was successful.
0357       *
0358       * \returns \c Success if computation was successful, \c NoConvergence otherwise.
0359       */
0360     EIGEN_DEVICE_FUNC
0361     ComputationInfo info() const
0362     {
0363       eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized.");
0364       return m_info;
0365     }
0366
0367     /** \brief Maximum number of iterations.
0368       *
0369       * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n
0370       * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK).
0371       */
0372     static const int m_maxIterations = 30;
0373
0374   protected:
0375     static EIGEN_DEVICE_FUNC
0376     void check_template_parameters()
0377     {
0378       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0379     }
0380
0381     EigenvectorsType m_eivec;
0382     RealVectorType m_eivalues;
0383     typename TridiagonalizationType::SubDiagonalType m_subdiag;
0384     typename TridiagonalizationType::CoeffVectorType m_hcoeffs;
0385     ComputationInfo m_info;
0386     bool m_isInitialized;
0387     bool m_eigenvectorsOk;
0388 };
0389
0390 namespace internal {
0391 /** \internal
0392   *
0393   * \eigenvalues_module \ingroup Eigenvalues_Module
0394   *
0395   * Performs a QR step on a tridiagonal symmetric matrix represented as a
0396   * pair of two vectors \a diag and \a subdiag.
0397   *
0398   * \param diag the diagonal part of the input selfadjoint tridiagonal matrix
0399   * \param subdiag the sub-diagonal part of the input selfadjoint tridiagonal matrix
0400   * \param start starting index of the submatrix to work on
0401   * \param end last+1 index of the submatrix to work on
0402   * \param matrixQ pointer to the column-major matrix holding the eigenvectors, can be 0
0403   * \param n size of the input matrix
0404   *
0405   * For compilation efficiency reasons, this procedure does not use eigen expression
0406   * for its arguments.
0407   *
0408   * Implemented from Golub's "Matrix Computations", algorithm 8.3.2:
0409   * "implicit symmetric QR step with Wilkinson shift"
0410   */
0411 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
0412 EIGEN_DEVICE_FUNC
0413 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n);
0414 }
0415
0416 template<typename MatrixType>
0417 template<typename InputType>
0418 EIGEN_DEVICE_FUNC
0419 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
0420 ::compute(const EigenBase<InputType>& a_matrix, int options)
0421 {
0422   check_template_parameters();
0423
0424   const InputType &matrix(a_matrix.derived());
0425
0426   EIGEN_USING_STD(abs);
0427   eigen_assert(matrix.cols() == matrix.rows());
0428   eigen_assert((options&~(EigVecMask|GenEigMask))==0
0429           && (options&EigVecMask)!=EigVecMask
0430           && "invalid option parameter");
0431   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
0432   Index n = matrix.cols();
0433   m_eivalues.resize(n,1);
0434
0435   if(n==1)
0436   {
0437     m_eivec = matrix;
0438     m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0));
0439     if(computeEigenvectors)
0440       m_eivec.setOnes(n,n);
0441     m_info = Success;
0442     m_isInitialized = true;
0443     m_eigenvectorsOk = computeEigenvectors;
0444     return *this;
0445   }
0446
0447   // declare some aliases
0448   RealVectorType& diag = m_eivalues;
0449   EigenvectorsType& mat = m_eivec;
0450
0451   // map the matrix coefficients to [-1:1] to avoid over- and underflow.
0452   mat = matrix.template triangularView<Lower>();
0453   RealScalar scale = mat.cwiseAbs().maxCoeff();
0454   if(scale==RealScalar(0)) scale = RealScalar(1);
0455   mat.template triangularView<Lower>() /= scale;
0456   m_subdiag.resize(n-1);
0457   m_hcoeffs.resize(n-1);
0458   internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors);
0459
0460   m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
0461
0462   // scale back the eigen values
0463   m_eivalues *= scale;
0464
0465   m_isInitialized = true;
0466   m_eigenvectorsOk = computeEigenvectors;
0467   return *this;
0468 }
0469
0470 template<typename MatrixType>
0471 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
0472 ::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options)
0473 {
0474   //TODO : Add an option to scale the values beforehand
0475   bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
0476
0477   m_eivalues = diag;
0478   m_subdiag = subdiag;
0479   if (computeEigenvectors)
0480   {
0481     m_eivec.setIdentity(diag.size(), diag.size());
0482   }
0483   m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec);
0484
0485   m_isInitialized = true;
0486   m_eigenvectorsOk = computeEigenvectors;
0487   return *this;
0488 }
0489
0490 namespace internal {
0491 /**
0492   * \internal
0493   * \brief Compute the eigendecomposition from a tridiagonal matrix
0494   *
0495   * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues
0496   * \param[in,out] subdiag : The subdiagonal part of the matrix (entries are modified during the decomposition)
0497   * \param[in] maxIterations : the maximum number of iterations
0498   * \param[in] computeEigenvectors : whether the eigenvectors have to be computed or not
0499   * \param[out] eivec : The matrix to store the eigenvectors if computeEigenvectors==true. Must be allocated on input.
0500   * \returns \c Success or \c NoConvergence
0501   */
0502 template<typename MatrixType, typename DiagType, typename SubDiagType>
0503 EIGEN_DEVICE_FUNC
0504 ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec)
0505 {
0506   ComputationInfo info;
0507   typedef typename MatrixType::Scalar Scalar;
0508
0509   Index n = diag.size();
0510   Index end = n-1;
0511   Index start = 0;
0512   Index iter = 0; // total number of iterations
0513
0514   typedef typename DiagType::RealScalar RealScalar;
0515   const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)();
0516   const RealScalar precision_inv = RealScalar(1)/NumTraits<RealScalar>::epsilon();
0517   while (end>0)
0518   {
0519     for (Index i = start; i<end; ++i) {
0520       if (numext::abs(subdiag[i]) < considerAsZero) {
0521         subdiag[i] = RealScalar(0);
0522       } else {
0523         // abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1]))
0524         // Scaled to prevent underflows.
0525         const RealScalar scaled_subdiag = precision_inv * subdiag[i];
0526         if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i])+numext::abs(diag[i+1]))) {
0527           subdiag[i] = RealScalar(0);
0528         }
0529       }
0530     }
0531
0532     // find the largest unreduced block at the end of the matrix.
0533     while (end>0 && subdiag[end-1]==RealScalar(0))
0534     {
0535       end--;
0536     }
0537     if (end<=0)
0538       break;
0539
0540     // if we spent too many iterations, we give up
0541     iter++;
0542     if(iter > maxIterations * n) break;
0543
0544     start = end - 1;
0545     while (start>0 && subdiag[start-1]!=0)
0546       start--;
0547
0548     internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n);
0549   }
0550   if (iter <= maxIterations * n)
0551     info = Success;
0552   else
0553     info = NoConvergence;
0554
0555   // Sort eigenvalues and corresponding vectors.
0556   // TODO make the sort optional ?
0557   // TODO use a better sort algorithm !!
0558   if (info == Success)
0559   {
0560     for (Index i = 0; i < n-1; ++i)
0561     {
0562       Index k;
0563       diag.segment(i,n-i).minCoeff(&k);
0564       if (k > 0)
0565       {
0566         numext::swap(diag[i], diag[k+i]);
0567         if(computeEigenvectors)
0568           eivec.col(i).swap(eivec.col(k+i));
0569       }
0570     }
0571   }
0572   return info;
0573 }
0574
0575 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues
0576 {
0577   EIGEN_DEVICE_FUNC
0578   static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options)
0579   { eig.compute(A,options); }
0580 };
0581
0582 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false>
0583 {
0584   typedef typename SolverType::MatrixType MatrixType;
0585   typedef typename SolverType::RealVectorType VectorType;
0586   typedef typename SolverType::Scalar Scalar;
0587   typedef typename SolverType::EigenvectorsType EigenvectorsType;
0588
0589
0590   /** \internal
0591    * Computes the roots of the characteristic polynomial of \a m.
0592    * For numerical stability m.trace() should be near zero and to avoid over- or underflow m should be normalized.
0593    */
0594   EIGEN_DEVICE_FUNC
0595   static inline void computeRoots(const MatrixType& m, VectorType& roots)
0596   {
0597     EIGEN_USING_STD(sqrt)
0598     EIGEN_USING_STD(atan2)
0599     EIGEN_USING_STD(cos)
0600     EIGEN_USING_STD(sin)
0601     const Scalar s_inv3 = Scalar(1)/Scalar(3);
0602     const Scalar s_sqrt3 = sqrt(Scalar(3));
0603
0604     // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0.  The
0605     // eigenvalues are the roots to this equation, all guaranteed to be
0606     // real-valued, because the matrix is symmetric.
0607     Scalar c0 = m(0,0)*m(1,1)*m(2,2) + Scalar(2)*m(1,0)*m(2,0)*m(2,1) - m(0,0)*m(2,1)*m(2,1) - m(1,1)*m(2,0)*m(2,0) - m(2,2)*m(1,0)*m(1,0);
0608     Scalar c1 = m(0,0)*m(1,1) - m(1,0)*m(1,0) + m(0,0)*m(2,2) - m(2,0)*m(2,0) + m(1,1)*m(2,2) - m(2,1)*m(2,1);
0609     Scalar c2 = m(0,0) + m(1,1) + m(2,2);
0610
0611     // Construct the parameters used in classifying the roots of the equation
0612     // and in solving the equation for the roots in closed form.
0613     Scalar c2_over_3 = c2*s_inv3;
0614     Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3;
0615     a_over_3 = numext::maxi(a_over_3, Scalar(0));
0616
0617     Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1));
0618
0619     Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b;
0620     q = numext::maxi(q, Scalar(0));
0621
0622     // Compute the eigenvalues by solving for the roots of the polynomial.
0623     Scalar rho = sqrt(a_over_3);
0624     Scalar theta = atan2(sqrt(q),half_b)*s_inv3;  // since sqrt(q) > 0, atan2 is in [0, pi] and theta is in [0, pi/3]
0625     Scalar cos_theta = cos(theta);
0626     Scalar sin_theta = sin(theta);
0627     // roots are already sorted, since cos is monotonically decreasing on [0, pi]
0628     roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3)
0629     roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3)
0630     roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta;
0631   }
0632
0633   EIGEN_DEVICE_FUNC
0634   static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative)
0635   {
0636     EIGEN_USING_STD(abs);
0637     EIGEN_USING_STD(sqrt);
0638     Index i0;
0639     // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal):
0640     mat.diagonal().cwiseAbs().maxCoeff(&i0);
0641     // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector,
0642     // so let's save it:
0643     representative = mat.col(i0);
0644     Scalar n0, n1;
0645     VectorType c0, c1;
0646     n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm();
0647     n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm();
0648     if(n0>n1) res = c0/sqrt(n0);
0649     else      res = c1/sqrt(n1);
0650
0651     return true;
0652   }
0653
0654   EIGEN_DEVICE_FUNC
0655   static inline void run(SolverType& solver, const MatrixType& mat, int options)
0656   {
0657     eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows());
0658     eigen_assert((options&~(EigVecMask|GenEigMask))==0
0659             && (options&EigVecMask)!=EigVecMask
0660             && "invalid option parameter");
0661     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
0662
0663     EigenvectorsType& eivecs = solver.m_eivec;
0664     VectorType& eivals = solver.m_eivalues;
0665
0666     // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
0667     Scalar shift = mat.trace() / Scalar(3);
0668     // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later
0669     MatrixType scaledMat = mat.template selfadjointView<Lower>();
0670     scaledMat.diagonal().array() -= shift;
0671     Scalar scale = scaledMat.cwiseAbs().maxCoeff();
0672     if(scale > 0) scaledMat /= scale;   // TODO for scale==0 we could save the remaining operations
0673
0674     // compute the eigenvalues
0675     computeRoots(scaledMat,eivals);
0676
0677     // compute the eigenvectors
0678     if(computeEigenvectors)
0679     {
0680       if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon())
0681       {
0682         // All three eigenvalues are numerically the same
0683         eivecs.setIdentity();
0684       }
0685       else
0686       {
0687         MatrixType tmp;
0688         tmp = scaledMat;
0689
0690         // Compute the eigenvector of the most distinct eigenvalue
0691         Scalar d0 = eivals(2) - eivals(1);
0692         Scalar d1 = eivals(1) - eivals(0);
0693         Index k(0), l(2);
0694         if(d0 > d1)
0695         {
0696           numext::swap(k,l);
0697           d0 = d1;
0698         }
0699
0700         // Compute the eigenvector of index k
0701         {
0702           tmp.diagonal().array () -= eivals(k);
0703           // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector.
0704           extract_kernel(tmp, eivecs.col(k), eivecs.col(l));
0705         }
0706
0707         // Compute eigenvector of index l
0708         if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1)
0709         {
0710           // If d0 is too small, then the two other eigenvalues are numerically the same,
0711           // and thus we only have to ortho-normalize the near orthogonal vector we saved above.
0712           eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l);
0713           eivecs.col(l).normalize();
0714         }
0715         else
0716         {
0717           tmp = scaledMat;
0718           tmp.diagonal().array () -= eivals(l);
0719
0720           VectorType dummy;
0721           extract_kernel(tmp, eivecs.col(l), dummy);
0722         }
0723
0724         // Compute last eigenvector from the other two
0725         eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized();
0726       }
0727     }
0728
0729     // Rescale back to the original size.
0730     eivals *= scale;
0731     eivals.array() += shift;
0732
0733     solver.m_info = Success;
0734     solver.m_isInitialized = true;
0735     solver.m_eigenvectorsOk = computeEigenvectors;
0736   }
0737 };
0738
0739 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel
0740 template<typename SolverType>
0741 struct direct_selfadjoint_eigenvalues<SolverType,2,false>
0742 {
0743   typedef typename SolverType::MatrixType MatrixType;
0744   typedef typename SolverType::RealVectorType VectorType;
0745   typedef typename SolverType::Scalar Scalar;
0746   typedef typename SolverType::EigenvectorsType EigenvectorsType;
0747
0748   EIGEN_DEVICE_FUNC
0749   static inline void computeRoots(const MatrixType& m, VectorType& roots)
0750   {
0751     EIGEN_USING_STD(sqrt);
0752     const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0)));
0753     const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1));
0754     roots(0) = t1 - t0;
0755     roots(1) = t1 + t0;
0756   }
0757
0758   EIGEN_DEVICE_FUNC
0759   static inline void run(SolverType& solver, const MatrixType& mat, int options)
0760   {
0761     EIGEN_USING_STD(sqrt);
0762     EIGEN_USING_STD(abs);
0763
0764     eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows());
0765     eigen_assert((options&~(EigVecMask|GenEigMask))==0
0766             && (options&EigVecMask)!=EigVecMask
0767             && "invalid option parameter");
0768     bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors;
0769
0770     EigenvectorsType& eivecs = solver.m_eivec;
0771     VectorType& eivals = solver.m_eivalues;
0772
0773     // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow.
0774     Scalar shift = mat.trace() / Scalar(2);
0775     MatrixType scaledMat = mat;
0776     scaledMat.coeffRef(0,1) = mat.coeff(1,0);
0777     scaledMat.diagonal().array() -= shift;
0778     Scalar scale = scaledMat.cwiseAbs().maxCoeff();
0779     if(scale > Scalar(0))
0780       scaledMat /= scale;
0781
0782     // Compute the eigenvalues
0783     computeRoots(scaledMat,eivals);
0784
0785     // compute the eigen vectors
0786     if(computeEigenvectors)
0787     {
0788       if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon())
0789       {
0790         eivecs.setIdentity();
0791       }
0792       else
0793       {
0794         scaledMat.diagonal().array () -= eivals(1);
0795         Scalar a2 = numext::abs2(scaledMat(0,0));
0796         Scalar c2 = numext::abs2(scaledMat(1,1));
0797         Scalar b2 = numext::abs2(scaledMat(1,0));
0798         if(a2>c2)
0799         {
0800           eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0);
0801           eivecs.col(1) /= sqrt(a2+b2);
0802         }
0803         else
0804         {
0805           eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0);
0806           eivecs.col(1) /= sqrt(c2+b2);
0807         }
0808
0809         eivecs.col(0) << eivecs.col(1).unitOrthogonal();
0810       }
0811     }
0812
0813     // Rescale back to the original size.
0814     eivals *= scale;
0815     eivals.array() += shift;
0816
0817     solver.m_info = Success;
0818     solver.m_isInitialized = true;
0819     solver.m_eigenvectorsOk = computeEigenvectors;
0820   }
0821 };
0822
0823 }
0824
0825 template<typename MatrixType>
0826 EIGEN_DEVICE_FUNC
0827 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType>
0828 ::computeDirect(const MatrixType& matrix, int options)
0829 {
0830   internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options);
0831   return *this;
0832 }
0833
0834 namespace internal {
0835
0836 // Francis implicit QR step.
0837 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index>
0838 EIGEN_DEVICE_FUNC
0839 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n)
0840 {
0841   // Wilkinson Shift.
0842   RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5);
0843   RealScalar e = subdiag[end-1];
0844   // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still
0845   // underflow thus leading to inf/NaN values when using the following commented code:
0846   //   RealScalar e2 = numext::abs2(subdiag[end-1]);
0847   //   RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2));
0848   // This explain the following, somewhat more complicated, version:
0849   RealScalar mu = diag[end];
0850   if(td==RealScalar(0)) {
0851     mu -= numext::abs(e);
0852   } else if (e != RealScalar(0)) {
0853     const RealScalar e2 = numext::abs2(e);
0854     const RealScalar h = numext::hypot(td,e);
0855     if(e2 == RealScalar(0)) {
0856       mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e);
0857     } else {
0858       mu -= e2 / (td + (td>RealScalar(0) ? h : -h));
0859     }
0860   }
0861
0862   RealScalar x = diag[start] - mu;
0863   RealScalar z = subdiag[start];
0864   // If z ever becomes zero, the Givens rotation will be the identity and
0865   // z will stay zero for all future iterations.
0866   for (Index k = start; k < end && z != RealScalar(0); ++k)
0867   {
0868     JacobiRotation<RealScalar> rot;
0869     rot.makeGivens(x, z);
0870
0871     // do T = G' T G
0872     RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k];
0873     RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1];
0874
0875     diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]);
0876     diag[k+1] = rot.s() * sdk + rot.c() * dkp1;
0877     subdiag[k] = rot.c() * sdk - rot.s() * dkp1;
0878
0879     if (k > start)
0880       subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z;
0881
0882     // "Chasing the bulge" to return to triangular form.
0883     x = subdiag[k];
0884     if (k < end - 1)
0885     {
0886       z = -rot.s() * subdiag[k+1];
0887       subdiag[k + 1] = rot.c() * subdiag[k+1];
0888     }
0889
0890     // apply the givens rotation to the unit matrix Q = Q * G
0891     if (matrixQ)
0892     {
0893       // FIXME if StorageOrder == RowMajor this operation is not very efficient
0894       Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n);
0895       q.applyOnTheRight(k,k+1,rot);
0896     }
0897   }
0898 }
0899
0900 } // end namespace internal
0901
0902 } // end namespace Eigen
0903
0904 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H