| // This file is part of Eigen, a lightweight C++ template library |
| // for linear algebra. |
| // |
| // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> |
| // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> |
| // |
| // This Source Code Form is subject to the terms of the Mozilla |
| // Public License v. 2.0. If a copy of the MPL was not distributed |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. |
| |
| #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H |
| #define EIGEN_SELFADJOINTEIGENSOLVER_H |
| |
| #include "./Tridiagonalization.h" |
| |
| #include "./InternalHeaderCheck.h" |
| |
| namespace Eigen { |
| |
| template<typename MatrixType_> |
| class GeneralizedSelfAdjointEigenSolver; |
| |
| namespace internal { |
| template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues; |
| |
| template<typename MatrixType, typename DiagType, typename SubDiagType> |
| EIGEN_DEVICE_FUNC |
| ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec); |
| } |
| |
| /** \eigenvalues_module \ingroup Eigenvalues_Module |
| * |
| * |
| * \class SelfAdjointEigenSolver |
| * |
| * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices |
| * |
| * \tparam MatrixType_ the type of the matrix of which we are computing the |
| * eigendecomposition; this is expected to be an instantiation of the Matrix |
| * class template. |
| * |
| * A matrix \f$ A \f$ is selfadjoint if it equals its adjoint. For real |
| * matrices, this means that the matrix is symmetric: it equals its |
| * transpose. This class computes the eigenvalues and eigenvectors of a |
| * selfadjoint matrix. These are the scalars \f$ \lambda \f$ and vectors |
| * \f$ v \f$ such that \f$ Av = \lambda v \f$. The eigenvalues of a |
| * selfadjoint matrix are always real. If \f$ D \f$ is a diagonal matrix with |
| * the eigenvalues on the diagonal, and \f$ V \f$ is a matrix with the |
| * eigenvectors as its columns, then \f$ A = V D V^{-1} \f$. This is called the |
| * eigendecomposition. |
| * |
| * For a selfadjoint matrix, \f$ V \f$ is unitary, meaning its inverse is equal |
| * to its adjoint, \f$ V^{-1} = V^{\dagger} \f$. If \f$ A \f$ is real, then |
| * \f$ V \f$ is also real and therefore orthogonal, meaning its inverse is |
| * equal to its transpose, \f$ V^{-1} = V^T \f$. |
| * |
| * The algorithm exploits the fact that the matrix is selfadjoint, making it |
| * faster and more accurate than the general purpose eigenvalue algorithms |
| * implemented in EigenSolver and ComplexEigenSolver. |
| * |
| * Only the \b lower \b triangular \b part of the input matrix is referenced. |
| * |
| * Call the function compute() to compute the eigenvalues and eigenvectors of |
| * a given matrix. Alternatively, you can use the |
| * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes |
| * the eigenvalues and eigenvectors at construction time. Once the eigenvalue |
| * and eigenvectors are computed, they can be retrieved with the eigenvalues() |
| * and eigenvectors() functions. |
| * |
| * The documentation for SelfAdjointEigenSolver(const MatrixType&, int) |
| * contains an example of the typical use of this class. |
| * |
| * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and |
| * the likes, see the class GeneralizedSelfAdjointEigenSolver. |
| * |
| * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver |
| */ |
| template<typename MatrixType_> class SelfAdjointEigenSolver |
| { |
| public: |
| |
| typedef MatrixType_ MatrixType; |
| enum { |
| Size = MatrixType::RowsAtCompileTime, |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, |
| Options = MatrixType::Options, |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime |
| }; |
| |
| /** \brief Scalar type for matrices of type \p MatrixType_. */ |
| typedef typename MatrixType::Scalar Scalar; |
| typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 |
| |
| typedef Matrix<Scalar,Size,Size,ColMajor,MaxColsAtCompileTime,MaxColsAtCompileTime> EigenvectorsType; |
| |
| /** \brief Real scalar type for \p MatrixType_. |
| * |
| * This is just \c Scalar if #Scalar is real (e.g., \c float or |
| * \c double), and the type of the real part of \c Scalar if #Scalar is |
| * complex. |
| */ |
| typedef typename NumTraits<Scalar>::Real RealScalar; |
| |
| friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>; |
| |
| /** \brief Type for vector of eigenvalues as returned by eigenvalues(). |
| * |
| * This is a column vector with entries of type #RealScalar. |
| * The length of the vector is the size of \p MatrixType_. |
| */ |
| typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; |
| typedef Tridiagonalization<MatrixType> TridiagonalizationType; |
| typedef typename TridiagonalizationType::SubDiagonalType SubDiagonalType; |
| |
| /** \brief Default constructor for fixed-size matrices. |
| * |
| * The default constructor is useful in cases in which the user intends to |
| * perform decompositions via compute(). This constructor |
| * can only be used if \p MatrixType_ is a fixed-size matrix; use |
| * SelfAdjointEigenSolver(Index) for dynamic-size matrices. |
| * |
| * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out |
| */ |
| EIGEN_DEVICE_FUNC |
| SelfAdjointEigenSolver() |
| : m_eivec(), |
| m_eivalues(), |
| m_subdiag(), |
| m_hcoeffs(), |
| m_info(InvalidInput), |
| m_isInitialized(false), |
| m_eigenvectorsOk(false) |
| { } |
| |
| /** \brief Constructor, pre-allocates memory for dynamic-size matrices. |
| * |
| * \param [in] size Positive integer, size of the matrix whose |
| * eigenvalues and eigenvectors will be computed. |
| * |
| * This constructor is useful for dynamic-size matrices, when the user |
| * intends to perform decompositions via compute(). The \p size |
| * parameter is only used as a hint. It is not an error to give a wrong |
| * \p size, but it may impair performance. |
| * |
| * \sa compute() for an example |
| */ |
| EIGEN_DEVICE_FUNC |
| explicit SelfAdjointEigenSolver(Index size) |
| : m_eivec(size, size), |
| m_eivalues(size), |
| m_subdiag(size > 1 ? size - 1 : 1), |
| m_hcoeffs(size > 1 ? size - 1 : 1), |
| m_isInitialized(false), |
| m_eigenvectorsOk(false) |
| {} |
| |
| /** \brief Constructor; computes eigendecomposition of given matrix. |
| * |
| * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to |
| * be computed. Only the lower triangular part of the matrix is referenced. |
| * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. |
| * |
| * This constructor calls compute(const MatrixType&, int) to compute the |
| * eigenvalues of the matrix \p matrix. The eigenvectors are computed if |
| * \p options equals #ComputeEigenvectors. |
| * |
| * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out |
| * |
| * \sa compute(const MatrixType&, int) |
| */ |
| template<typename InputType> |
| EIGEN_DEVICE_FUNC |
| explicit SelfAdjointEigenSolver(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors) |
| : m_eivec(matrix.rows(), matrix.cols()), |
| m_eivalues(matrix.cols()), |
| m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), |
| m_hcoeffs(matrix.cols() > 1 ? matrix.cols() - 1 : 1), |
| m_isInitialized(false), |
| m_eigenvectorsOk(false) |
| { |
| compute(matrix.derived(), options); |
| } |
| |
| /** \brief Computes eigendecomposition of given matrix. |
| * |
| * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to |
| * be computed. Only the lower triangular part of the matrix is referenced. |
| * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. |
| * \returns Reference to \c *this |
| * |
| * This function computes the eigenvalues of \p matrix. The eigenvalues() |
| * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, |
| * then the eigenvectors are also computed and can be retrieved by |
| * calling eigenvectors(). |
| * |
| * This implementation uses a symmetric QR algorithm. The matrix is first |
| * reduced to tridiagonal form using the Tridiagonalization class. The |
| * tridiagonal matrix is then brought to diagonal form with implicit |
| * symmetric QR steps with Wilkinson shift. Details can be found in |
| * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>. |
| * |
| * The cost of the computation is about \f$ 9n^3 \f$ if the eigenvectors |
| * are required and \f$ 4n^3/3 \f$ if they are not required. |
| * |
| * This method reuses the memory in the SelfAdjointEigenSolver object that |
| * was allocated when the object was constructed, if the size of the |
| * matrix does not change. |
| * |
| * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out |
| * |
| * \sa SelfAdjointEigenSolver(const MatrixType&, int) |
| */ |
| template<typename InputType> |
| EIGEN_DEVICE_FUNC |
| SelfAdjointEigenSolver& compute(const EigenBase<InputType>& matrix, int options = ComputeEigenvectors); |
| |
| /** \brief Computes eigendecomposition of given matrix using a closed-form algorithm |
| * |
| * This is a variant of compute(const MatrixType&, int options) which |
| * directly solves the underlying polynomial equation. |
| * |
| * Currently only 2x2 and 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). |
| * |
| * This method is usually significantly faster than the QR iterative algorithm |
| * but it might also be less accurate. It is also worth noting that |
| * for 3x3 matrices it involves trigonometric operations which are |
| * not necessarily available for all scalar types. |
| * |
| * For the 3x3 case, we observed the following worst case relative error regarding the eigenvalues: |
| * - double: 1e-8 |
| * - float: 1e-3 |
| * |
| * \sa compute(const MatrixType&, int options) |
| */ |
| EIGEN_DEVICE_FUNC |
| SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); |
| |
| /** |
| *\brief Computes the eigen decomposition from a tridiagonal symmetric matrix |
| * |
| * \param[in] diag The vector containing the diagonal of the matrix. |
| * \param[in] subdiag The subdiagonal of the matrix. |
| * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. |
| * \returns Reference to \c *this |
| * |
| * This function assumes that the matrix has been reduced to tridiagonal form. |
| * |
| * \sa compute(const MatrixType&, int) for more information |
| */ |
| SelfAdjointEigenSolver& computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options=ComputeEigenvectors); |
| |
| /** \brief Returns the eigenvectors of given matrix. |
| * |
| * \returns A const reference to the matrix whose columns are the eigenvectors. |
| * |
| * \pre The eigenvectors have been computed before. |
| * |
| * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding |
| * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The |
| * eigenvectors are normalized to have (Euclidean) norm equal to one. If |
| * this object was used to solve the eigenproblem for the selfadjoint |
| * matrix \f$ A \f$, then the matrix returned by this function is the |
| * matrix \f$ V \f$ in the eigendecomposition \f$ A = V D V^{-1} \f$. |
| * |
| * For a selfadjoint matrix, \f$ V \f$ is unitary, meaning its inverse is equal |
| * to its adjoint, \f$ V^{-1} = V^{\dagger} \f$. If \f$ A \f$ is real, then |
| * \f$ V \f$ is also real and therefore orthogonal, meaning its inverse is |
| * equal to its transpose, \f$ V^{-1} = V^T \f$. |
| * |
| * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out |
| * |
| * \sa eigenvalues() |
| */ |
| EIGEN_DEVICE_FUNC |
| const EigenvectorsType& eigenvectors() const |
| { |
| eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); |
| eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
| return m_eivec; |
| } |
| |
| /** \brief Returns the eigenvalues of given matrix. |
| * |
| * \returns A const reference to the column vector containing the eigenvalues. |
| * |
| * \pre The eigenvalues have been computed before. |
| * |
| * The eigenvalues are repeated according to their algebraic multiplicity, |
| * so there are as many eigenvalues as rows in the matrix. The eigenvalues |
| * are sorted in increasing order. |
| * |
| * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out |
| * |
| * \sa eigenvectors(), MatrixBase::eigenvalues() |
| */ |
| EIGEN_DEVICE_FUNC |
| const RealVectorType& eigenvalues() const |
| { |
| eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); |
| return m_eivalues; |
| } |
| |
| /** \brief Computes the positive-definite square root of the matrix. |
| * |
| * \returns the positive-definite square root of the matrix |
| * |
| * \pre The eigenvalues and eigenvectors of a positive-definite matrix |
| * have been computed before. |
| * |
| * The square root of a positive-definite matrix \f$ A \f$ is the |
| * positive-definite matrix whose square equals \f$ A \f$. This function |
| * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the |
| * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$. |
| * |
| * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out |
| * |
| * \sa operatorInverseSqrt(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> |
| */ |
| EIGEN_DEVICE_FUNC |
| MatrixType operatorSqrt() const |
| { |
| eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); |
| eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
| return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); |
| } |
| |
| /** \brief Computes the inverse square root of the matrix. |
| * |
| * \returns the inverse positive-definite square root of the matrix |
| * |
| * \pre The eigenvalues and eigenvectors of a positive-definite matrix |
| * have been computed before. |
| * |
| * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to |
| * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is |
| * cheaper than first computing the square root with operatorSqrt() and |
| * then its inverse with MatrixBase::inverse(). |
| * |
| * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp |
| * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out |
| * |
| * \sa operatorSqrt(), MatrixBase::inverse(), <a href="unsupported/group__MatrixFunctions__Module.html">MatrixFunctions Module</a> |
| */ |
| EIGEN_DEVICE_FUNC |
| MatrixType operatorInverseSqrt() const |
| { |
| eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); |
| eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); |
| return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); |
| } |
| |
| /** \brief Reports whether previous computation was successful. |
| * |
| * \returns \c Success if computation was successful, \c NoConvergence otherwise. |
| */ |
| EIGEN_DEVICE_FUNC |
| ComputationInfo info() const |
| { |
| eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); |
| return m_info; |
| } |
| |
| /** \brief Maximum number of iterations. |
| * |
| * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n |
| * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK). |
| */ |
| static const int m_maxIterations = 30; |
| |
| protected: |
| EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar) |
| |
| EigenvectorsType m_eivec; |
| RealVectorType m_eivalues; |
| typename TridiagonalizationType::SubDiagonalType m_subdiag; |
| typename TridiagonalizationType::CoeffVectorType m_hcoeffs; |
| ComputationInfo m_info; |
| bool m_isInitialized; |
| bool m_eigenvectorsOk; |
| }; |
| |
| namespace internal { |
| /** \internal |
| * |
| * \eigenvalues_module \ingroup Eigenvalues_Module |
| * |
| * Performs a QR step on a tridiagonal symmetric matrix represented as a |
| * pair of two vectors \a diag and \a subdiag. |
| * |
| * \param diag the diagonal part of the input selfadjoint tridiagonal matrix |
| * \param subdiag the sub-diagonal part of the input selfadjoint tridiagonal matrix |
| * \param start starting index of the submatrix to work on |
| * \param end last+1 index of the submatrix to work on |
| * \param matrixQ pointer to the column-major matrix holding the eigenvectors, can be 0 |
| * \param n size of the input matrix |
| * |
| * For compilation efficiency reasons, this procedure does not use eigen expression |
| * for its arguments. |
| * |
| * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: |
| * "implicit symmetric QR step with Wilkinson shift" |
| */ |
| template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> |
| EIGEN_DEVICE_FUNC |
| static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); |
| } |
| |
| template<typename MatrixType> |
| template<typename InputType> |
| EIGEN_DEVICE_FUNC |
| SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> |
| ::compute(const EigenBase<InputType>& a_matrix, int options) |
| { |
| const InputType &matrix(a_matrix.derived()); |
| |
| EIGEN_USING_STD(abs); |
| eigen_assert(matrix.cols() == matrix.rows()); |
| eigen_assert((options&~(EigVecMask|GenEigMask))==0 |
| && (options&EigVecMask)!=EigVecMask |
| && "invalid option parameter"); |
| bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; |
| Index n = matrix.cols(); |
| m_eivalues.resize(n,1); |
| |
| if(n==1) |
| { |
| m_eivec = matrix; |
| m_eivalues.coeffRef(0,0) = numext::real(m_eivec.coeff(0,0)); |
| if(computeEigenvectors) |
| m_eivec.setOnes(n,n); |
| m_info = Success; |
| m_isInitialized = true; |
| m_eigenvectorsOk = computeEigenvectors; |
| return *this; |
| } |
| |
| // declare some aliases |
| RealVectorType& diag = m_eivalues; |
| EigenvectorsType& mat = m_eivec; |
| |
| // map the matrix coefficients to [-1:1] to avoid over- and underflow. |
| mat = matrix.template triangularView<Lower>(); |
| RealScalar scale = mat.cwiseAbs().maxCoeff(); |
| if(numext::is_exactly_zero(scale)) scale = RealScalar(1); |
| mat.template triangularView<Lower>() /= scale; |
| m_subdiag.resize(n-1); |
| m_hcoeffs.resize(n-1); |
| internal::tridiagonalization_inplace(mat, diag, m_subdiag, m_hcoeffs, computeEigenvectors); |
| |
| m_info = internal::computeFromTridiagonal_impl(diag, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); |
| |
| // scale back the eigen values |
| m_eivalues *= scale; |
| |
| m_isInitialized = true; |
| m_eigenvectorsOk = computeEigenvectors; |
| return *this; |
| } |
| |
| template<typename MatrixType> |
| SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> |
| ::computeFromTridiagonal(const RealVectorType& diag, const SubDiagonalType& subdiag , int options) |
| { |
| //TODO : Add an option to scale the values beforehand |
| bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; |
| |
| m_eivalues = diag; |
| m_subdiag = subdiag; |
| if (computeEigenvectors) |
| { |
| m_eivec.setIdentity(diag.size(), diag.size()); |
| } |
| m_info = internal::computeFromTridiagonal_impl(m_eivalues, m_subdiag, m_maxIterations, computeEigenvectors, m_eivec); |
| |
| m_isInitialized = true; |
| m_eigenvectorsOk = computeEigenvectors; |
| return *this; |
| } |
| |
| namespace internal { |
| /** |
| * \internal |
| * \brief Compute the eigendecomposition from a tridiagonal matrix |
| * |
| * \param[in,out] diag : On input, the diagonal of the matrix, on output the eigenvalues |
| * \param[in,out] subdiag : The subdiagonal part of the matrix (entries are modified during the decomposition) |
| * \param[in] maxIterations : the maximum number of iterations |
| * \param[in] computeEigenvectors : whether the eigenvectors have to be computed or not |
| * \param[out] eivec : The matrix to store the eigenvectors if computeEigenvectors==true. Must be allocated on input. |
| * \returns \c Success or \c NoConvergence |
| */ |
| template<typename MatrixType, typename DiagType, typename SubDiagType> |
| EIGEN_DEVICE_FUNC |
| ComputationInfo computeFromTridiagonal_impl(DiagType& diag, SubDiagType& subdiag, const Index maxIterations, bool computeEigenvectors, MatrixType& eivec) |
| { |
| ComputationInfo info; |
| typedef typename MatrixType::Scalar Scalar; |
| |
| Index n = diag.size(); |
| Index end = n-1; |
| Index start = 0; |
| Index iter = 0; // total number of iterations |
| |
| typedef typename DiagType::RealScalar RealScalar; |
| const RealScalar considerAsZero = (std::numeric_limits<RealScalar>::min)(); |
| const RealScalar precision_inv = RealScalar(1)/NumTraits<RealScalar>::epsilon(); |
| while (end>0) |
| { |
| for (Index i = start; i<end; ++i) { |
| if (numext::abs(subdiag[i]) < considerAsZero) { |
| subdiag[i] = RealScalar(0); |
| } else { |
| // abs(subdiag[i]) <= epsilon * sqrt(abs(diag[i]) + abs(diag[i+1])) |
| // Scaled to prevent underflows. |
| const RealScalar scaled_subdiag = precision_inv * subdiag[i]; |
| if (scaled_subdiag * scaled_subdiag <= (numext::abs(diag[i])+numext::abs(diag[i+1]))) { |
| subdiag[i] = RealScalar(0); |
| } |
| } |
| } |
| |
| // find the largest unreduced block at the end of the matrix. |
| while (end>0 && numext::is_exactly_zero(subdiag[end - 1])) |
| { |
| end--; |
| } |
| if (end<=0) |
| break; |
| |
| // if we spent too many iterations, we give up |
| iter++; |
| if(iter > maxIterations * n) break; |
| |
| start = end - 1; |
| while (start>0 && !numext::is_exactly_zero(subdiag[start - 1])) |
| start--; |
| |
| internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), subdiag.data(), start, end, computeEigenvectors ? eivec.data() : (Scalar*)0, n); |
| } |
| if (iter <= maxIterations * n) |
| info = Success; |
| else |
| info = NoConvergence; |
| |
| // Sort eigenvalues and corresponding vectors. |
| // TODO make the sort optional ? |
| // TODO use a better sort algorithm !! |
| if (info == Success) |
| { |
| for (Index i = 0; i < n-1; ++i) |
| { |
| Index k; |
| diag.segment(i,n-i).minCoeff(&k); |
| if (k > 0) |
| { |
| numext::swap(diag[i], diag[k+i]); |
| if(computeEigenvectors) |
| eivec.col(i).swap(eivec.col(k+i)); |
| } |
| } |
| } |
| return info; |
| } |
| |
| template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues |
| { |
| EIGEN_DEVICE_FUNC |
| static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) |
| { eig.compute(A,options); } |
| }; |
| |
| template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false> |
| { |
| typedef typename SolverType::MatrixType MatrixType; |
| typedef typename SolverType::RealVectorType VectorType; |
| typedef typename SolverType::Scalar Scalar; |
| typedef typename SolverType::EigenvectorsType EigenvectorsType; |
| |
| |
| /** \internal |
| * Computes the roots of the characteristic polynomial of \a m. |
| * For numerical stability m.trace() should be near zero and to avoid over- or underflow m should be normalized. |
| */ |
| EIGEN_DEVICE_FUNC |
| static inline void computeRoots(const MatrixType& m, VectorType& roots) |
| { |
| EIGEN_USING_STD(sqrt) |
| EIGEN_USING_STD(atan2) |
| EIGEN_USING_STD(cos) |
| EIGEN_USING_STD(sin) |
| const Scalar s_inv3 = Scalar(1)/Scalar(3); |
| const Scalar s_sqrt3 = sqrt(Scalar(3)); |
| |
| // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The |
| // eigenvalues are the roots to this equation, all guaranteed to be |
| // real-valued, because the matrix is symmetric. |
| 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); |
| 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); |
| Scalar c2 = m(0,0) + m(1,1) + m(2,2); |
| |
| // Construct the parameters used in classifying the roots of the equation |
| // and in solving the equation for the roots in closed form. |
| Scalar c2_over_3 = c2*s_inv3; |
| Scalar a_over_3 = (c2*c2_over_3 - c1)*s_inv3; |
| a_over_3 = numext::maxi(a_over_3, Scalar(0)); |
| |
| Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); |
| |
| Scalar q = a_over_3*a_over_3*a_over_3 - half_b*half_b; |
| q = numext::maxi(q, Scalar(0)); |
| |
| // Compute the eigenvalues by solving for the roots of the polynomial. |
| Scalar rho = sqrt(a_over_3); |
| 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] |
| Scalar cos_theta = cos(theta); |
| Scalar sin_theta = sin(theta); |
| // roots are already sorted, since cos is monotonically decreasing on [0, pi] |
| roots(0) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); // == 2*rho*cos(theta+2pi/3) |
| roots(1) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); // == 2*rho*cos(theta+ pi/3) |
| roots(2) = c2_over_3 + Scalar(2)*rho*cos_theta; |
| } |
| |
| EIGEN_DEVICE_FUNC |
| static inline bool extract_kernel(MatrixType& mat, Ref<VectorType> res, Ref<VectorType> representative) |
| { |
| EIGEN_USING_STD(abs); |
| EIGEN_USING_STD(sqrt); |
| Index i0; |
| // Find non-zero column i0 (by construction, there must exist a non zero coefficient on the diagonal): |
| mat.diagonal().cwiseAbs().maxCoeff(&i0); |
| // mat.col(i0) is a good candidate for an orthogonal vector to the current eigenvector, |
| // so let's save it: |
| representative = mat.col(i0); |
| Scalar n0, n1; |
| VectorType c0, c1; |
| n0 = (c0 = representative.cross(mat.col((i0+1)%3))).squaredNorm(); |
| n1 = (c1 = representative.cross(mat.col((i0+2)%3))).squaredNorm(); |
| if(n0>n1) res = c0/sqrt(n0); |
| else res = c1/sqrt(n1); |
| |
| return true; |
| } |
| |
| EIGEN_DEVICE_FUNC |
| static inline void run(SolverType& solver, const MatrixType& mat, int options) |
| { |
| eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); |
| eigen_assert((options&~(EigVecMask|GenEigMask))==0 |
| && (options&EigVecMask)!=EigVecMask |
| && "invalid option parameter"); |
| bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; |
| |
| EigenvectorsType& eivecs = solver.m_eivec; |
| VectorType& eivals = solver.m_eivalues; |
| |
| // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. |
| Scalar shift = mat.trace() / Scalar(3); |
| // TODO Avoid this copy. Currently it is necessary to suppress bogus values when determining maxCoeff and for computing the eigenvectors later |
| MatrixType scaledMat = mat.template selfadjointView<Lower>(); |
| scaledMat.diagonal().array() -= shift; |
| Scalar scale = scaledMat.cwiseAbs().maxCoeff(); |
| if(scale > 0) scaledMat /= scale; // TODO for scale==0 we could save the remaining operations |
| |
| // compute the eigenvalues |
| computeRoots(scaledMat,eivals); |
| |
| // compute the eigenvectors |
| if(computeEigenvectors) |
| { |
| if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon()) |
| { |
| // All three eigenvalues are numerically the same |
| eivecs.setIdentity(); |
| } |
| else |
| { |
| MatrixType tmp; |
| tmp = scaledMat; |
| |
| // Compute the eigenvector of the most distinct eigenvalue |
| Scalar d0 = eivals(2) - eivals(1); |
| Scalar d1 = eivals(1) - eivals(0); |
| Index k(0), l(2); |
| if(d0 > d1) |
| { |
| numext::swap(k,l); |
| d0 = d1; |
| } |
| |
| // Compute the eigenvector of index k |
| { |
| tmp.diagonal().array () -= eivals(k); |
| // By construction, 'tmp' is of rank 2, and its kernel corresponds to the respective eigenvector. |
| extract_kernel(tmp, eivecs.col(k), eivecs.col(l)); |
| } |
| |
| // Compute eigenvector of index l |
| if(d0<=2*Eigen::NumTraits<Scalar>::epsilon()*d1) |
| { |
| // If d0 is too small, then the two other eigenvalues are numerically the same, |
| // and thus we only have to ortho-normalize the near orthogonal vector we saved above. |
| eivecs.col(l) -= eivecs.col(k).dot(eivecs.col(l))*eivecs.col(l); |
| eivecs.col(l).normalize(); |
| } |
| else |
| { |
| tmp = scaledMat; |
| tmp.diagonal().array () -= eivals(l); |
| |
| VectorType dummy; |
| extract_kernel(tmp, eivecs.col(l), dummy); |
| } |
| |
| // Compute last eigenvector from the other two |
| eivecs.col(1) = eivecs.col(2).cross(eivecs.col(0)).normalized(); |
| } |
| } |
| |
| // Rescale back to the original size. |
| eivals *= scale; |
| eivals.array() += shift; |
| |
| solver.m_info = Success; |
| solver.m_isInitialized = true; |
| solver.m_eigenvectorsOk = computeEigenvectors; |
| } |
| }; |
| |
| // 2x2 direct eigenvalues decomposition, code from Hauke Heibel |
| template<typename SolverType> |
| struct direct_selfadjoint_eigenvalues<SolverType,2,false> |
| { |
| typedef typename SolverType::MatrixType MatrixType; |
| typedef typename SolverType::RealVectorType VectorType; |
| typedef typename SolverType::Scalar Scalar; |
| typedef typename SolverType::EigenvectorsType EigenvectorsType; |
| |
| EIGEN_DEVICE_FUNC |
| static inline void computeRoots(const MatrixType& m, VectorType& roots) |
| { |
| EIGEN_USING_STD(sqrt); |
| const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*numext::abs2(m(1,0))); |
| const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1)); |
| roots(0) = t1 - t0; |
| roots(1) = t1 + t0; |
| } |
| |
| EIGEN_DEVICE_FUNC |
| static inline void run(SolverType& solver, const MatrixType& mat, int options) |
| { |
| EIGEN_USING_STD(sqrt); |
| EIGEN_USING_STD(abs); |
| |
| eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); |
| eigen_assert((options&~(EigVecMask|GenEigMask))==0 |
| && (options&EigVecMask)!=EigVecMask |
| && "invalid option parameter"); |
| bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; |
| |
| EigenvectorsType& eivecs = solver.m_eivec; |
| VectorType& eivals = solver.m_eivalues; |
| |
| // Shift the matrix to the mean eigenvalue and map the matrix coefficients to [-1:1] to avoid over- and underflow. |
| Scalar shift = mat.trace() / Scalar(2); |
| MatrixType scaledMat = mat; |
| scaledMat.coeffRef(0,1) = mat.coeff(1,0); |
| scaledMat.diagonal().array() -= shift; |
| Scalar scale = scaledMat.cwiseAbs().maxCoeff(); |
| if(scale > Scalar(0)) |
| scaledMat /= scale; |
| |
| // Compute the eigenvalues |
| computeRoots(scaledMat,eivals); |
| |
| // compute the eigen vectors |
| if(computeEigenvectors) |
| { |
| if((eivals(1)-eivals(0))<=abs(eivals(1))*Eigen::NumTraits<Scalar>::epsilon()) |
| { |
| eivecs.setIdentity(); |
| } |
| else |
| { |
| scaledMat.diagonal().array () -= eivals(1); |
| Scalar a2 = numext::abs2(scaledMat(0,0)); |
| Scalar c2 = numext::abs2(scaledMat(1,1)); |
| Scalar b2 = numext::abs2(scaledMat(1,0)); |
| if(a2>c2) |
| { |
| eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); |
| eivecs.col(1) /= sqrt(a2+b2); |
| } |
| else |
| { |
| eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); |
| eivecs.col(1) /= sqrt(c2+b2); |
| } |
| |
| eivecs.col(0) << eivecs.col(1).unitOrthogonal(); |
| } |
| } |
| |
| // Rescale back to the original size. |
| eivals *= scale; |
| eivals.array() += shift; |
| |
| solver.m_info = Success; |
| solver.m_isInitialized = true; |
| solver.m_eigenvectorsOk = computeEigenvectors; |
| } |
| }; |
| |
| } |
| |
| template<typename MatrixType> |
| EIGEN_DEVICE_FUNC |
| SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> |
| ::computeDirect(const MatrixType& matrix, int options) |
| { |
| internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options); |
| return *this; |
| } |
| |
| namespace internal { |
| |
| // Francis implicit QR step. |
| template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> |
| EIGEN_DEVICE_FUNC |
| static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) |
| { |
| // Wilkinson Shift. |
| RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); |
| RealScalar e = subdiag[end-1]; |
| // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still |
| // underflow thus leading to inf/NaN values when using the following commented code: |
| // RealScalar e2 = numext::abs2(subdiag[end-1]); |
| // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); |
| // This explain the following, somewhat more complicated, version: |
| RealScalar mu = diag[end]; |
| if(numext::is_exactly_zero(td)) { |
| mu -= numext::abs(e); |
| } else if (!numext::is_exactly_zero(e)) { |
| const RealScalar e2 = numext::abs2(e); |
| const RealScalar h = numext::hypot(td,e); |
| if(numext::is_exactly_zero(e2)) { |
| mu -= e / ((td + (td>RealScalar(0) ? h : -h)) / e); |
| } else { |
| mu -= e2 / (td + (td>RealScalar(0) ? h : -h)); |
| } |
| } |
| |
| RealScalar x = diag[start] - mu; |
| RealScalar z = subdiag[start]; |
| // If z ever becomes zero, the Givens rotation will be the identity and |
| // z will stay zero for all future iterations. |
| for (Index k = start; k < end && !numext::is_exactly_zero(z); ++k) |
| { |
| JacobiRotation<RealScalar> rot; |
| rot.makeGivens(x, z); |
| |
| // do T = G' T G |
| RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; |
| RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; |
| |
| diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); |
| diag[k+1] = rot.s() * sdk + rot.c() * dkp1; |
| subdiag[k] = rot.c() * sdk - rot.s() * dkp1; |
| |
| if (k > start) |
| subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; |
| |
| // "Chasing the bulge" to return to triangular form. |
| x = subdiag[k]; |
| if (k < end - 1) |
| { |
| z = -rot.s() * subdiag[k+1]; |
| subdiag[k + 1] = rot.c() * subdiag[k+1]; |
| } |
| |
| // apply the givens rotation to the unit matrix Q = Q * G |
| if (matrixQ) |
| { |
| // FIXME if StorageOrder == RowMajor this operation is not very efficient |
| Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n); |
| q.applyOnTheRight(k,k+1,rot); |
| } |
| } |
| } |
| |
| } // end namespace internal |
| |
| } // end namespace Eigen |
| |
| #endif // EIGEN_SELFADJOINTEIGENSOLVER_H |
