1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2010 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_SELFADJOINTEIGENSOLVER_H 12 #define EIGEN_SELFADJOINTEIGENSOLVER_H 13 14 #include "./Tridiagonalization.h" 15 16 namespace Eigen { 17 18 template<typename _MatrixType> 19 class GeneralizedSelfAdjointEigenSolver; 20 21 namespace internal { 22 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues; 23 } 24 25 /** \eigenvalues_module \ingroup Eigenvalues_Module 26 * 27 * 28 * \class SelfAdjointEigenSolver 29 * 30 * \brief Computes eigenvalues and eigenvectors of selfadjoint matrices 31 * 32 * \tparam _MatrixType the type of the matrix of which we are computing the 33 * eigendecomposition; this is expected to be an instantiation of the Matrix 34 * class template. 35 * 36 * A matrix \f$ A \f$ is selfadjoint if it equals its adjoint. For real 37 * matrices, this means that the matrix is symmetric: it equals its 38 * transpose. This class computes the eigenvalues and eigenvectors of a 39 * selfadjoint matrix. These are the scalars \f$ \lambda \f$ and vectors 40 * \f$ v \f$ such that \f$ Av = \lambda v \f$. The eigenvalues of a 41 * selfadjoint matrix are always real. If \f$ D \f$ is a diagonal matrix with 42 * the eigenvalues on the diagonal, and \f$ V \f$ is a matrix with the 43 * eigenvectors as its columns, then \f$ A = V D V^{-1} \f$ (for selfadjoint 44 * matrices, the matrix \f$ V \f$ is always invertible). This is called the 45 * eigendecomposition. 46 * 47 * The algorithm exploits the fact that the matrix is selfadjoint, making it 48 * faster and more accurate than the general purpose eigenvalue algorithms 49 * implemented in EigenSolver and ComplexEigenSolver. 50 * 51 * Only the \b lower \b triangular \b part of the input matrix is referenced. 52 * 53 * Call the function compute() to compute the eigenvalues and eigenvectors of 54 * a given matrix. Alternatively, you can use the 55 * SelfAdjointEigenSolver(const MatrixType&, int) constructor which computes 56 * the eigenvalues and eigenvectors at construction time. Once the eigenvalue 57 * and eigenvectors are computed, they can be retrieved with the eigenvalues() 58 * and eigenvectors() functions. 59 * 60 * The documentation for SelfAdjointEigenSolver(const MatrixType&, int) 61 * contains an example of the typical use of this class. 62 * 63 * To solve the \em generalized eigenvalue problem \f$ Av = \lambda Bv \f$ and 64 * the likes, see the class GeneralizedSelfAdjointEigenSolver. 65 * 66 * \sa MatrixBase::eigenvalues(), class EigenSolver, class ComplexEigenSolver 67 */ 68 template<typename _MatrixType> class SelfAdjointEigenSolver 69 { 70 public: 71 72 typedef _MatrixType MatrixType; 73 enum { 74 Size = MatrixType::RowsAtCompileTime, 75 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 76 Options = MatrixType::Options, 77 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 78 }; 79 80 /** \brief Scalar type for matrices of type \p _MatrixType. */ 81 typedef typename MatrixType::Scalar Scalar; 82 typedef typename MatrixType::Index Index; 83 84 /** \brief Real scalar type for \p _MatrixType. 85 * 86 * This is just \c Scalar if #Scalar is real (e.g., \c float or 87 * \c double), and the type of the real part of \c Scalar if #Scalar is 88 * complex. 89 */ 90 typedef typename NumTraits<Scalar>::Real RealScalar; 91 92 friend struct internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>; 93 94 /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 95 * 96 * This is a column vector with entries of type #RealScalar. 97 * The length of the vector is the size of \p _MatrixType. 98 */ 99 typedef typename internal::plain_col_type<MatrixType, RealScalar>::type RealVectorType; 100 typedef Tridiagonalization<MatrixType> TridiagonalizationType; 101 102 /** \brief Default constructor for fixed-size matrices. 103 * 104 * The default constructor is useful in cases in which the user intends to 105 * perform decompositions via compute(). This constructor 106 * can only be used if \p _MatrixType is a fixed-size matrix; use 107 * SelfAdjointEigenSolver(Index) for dynamic-size matrices. 108 * 109 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver.cpp 110 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver.out 111 */ 112 SelfAdjointEigenSolver() 113 : m_eivec(), 114 m_eivalues(), 115 m_subdiag(), 116 m_isInitialized(false) 117 { } 118 119 /** \brief Constructor, pre-allocates memory for dynamic-size matrices. 120 * 121 * \param [in] size Positive integer, size of the matrix whose 122 * eigenvalues and eigenvectors will be computed. 123 * 124 * This constructor is useful for dynamic-size matrices, when the user 125 * intends to perform decompositions via compute(). The \p size 126 * parameter is only used as a hint. It is not an error to give a wrong 127 * \p size, but it may impair performance. 128 * 129 * \sa compute() for an example 130 */ 131 SelfAdjointEigenSolver(Index size) 132 : m_eivec(size, size), 133 m_eivalues(size), 134 m_subdiag(size > 1 ? size - 1 : 1), 135 m_isInitialized(false) 136 {} 137 138 /** \brief Constructor; computes eigendecomposition of given matrix. 139 * 140 * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to 141 * be computed. Only the lower triangular part of the matrix is referenced. 142 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 143 * 144 * This constructor calls compute(const MatrixType&, int) to compute the 145 * eigenvalues of the matrix \p matrix. The eigenvectors are computed if 146 * \p options equals #ComputeEigenvectors. 147 * 148 * Example: \include SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.cpp 149 * Output: \verbinclude SelfAdjointEigenSolver_SelfAdjointEigenSolver_MatrixType.out 150 * 151 * \sa compute(const MatrixType&, int) 152 */ 153 SelfAdjointEigenSolver(const MatrixType& matrix, int options = ComputeEigenvectors) 154 : m_eivec(matrix.rows(), matrix.cols()), 155 m_eivalues(matrix.cols()), 156 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), 157 m_isInitialized(false) 158 { 159 compute(matrix, options); 160 } 161 162 /** \brief Computes eigendecomposition of given matrix. 163 * 164 * \param[in] matrix Selfadjoint matrix whose eigendecomposition is to 165 * be computed. Only the lower triangular part of the matrix is referenced. 166 * \param[in] options Can be #ComputeEigenvectors (default) or #EigenvaluesOnly. 167 * \returns Reference to \c *this 168 * 169 * This function computes the eigenvalues of \p matrix. The eigenvalues() 170 * function can be used to retrieve them. If \p options equals #ComputeEigenvectors, 171 * then the eigenvectors are also computed and can be retrieved by 172 * calling eigenvectors(). 173 * 174 * This implementation uses a symmetric QR algorithm. The matrix is first 175 * reduced to tridiagonal form using the Tridiagonalization class. The 176 * tridiagonal matrix is then brought to diagonal form with implicit 177 * symmetric QR steps with Wilkinson shift. Details can be found in 178 * Section 8.3 of Golub \& Van Loan, <i>%Matrix Computations</i>. 179 * 180 * The cost of the computation is about \f$ 9n^3 \f$ if the eigenvectors 181 * are required and \f$ 4n^3/3 \f$ if they are not required. 182 * 183 * This method reuses the memory in the SelfAdjointEigenSolver object that 184 * was allocated when the object was constructed, if the size of the 185 * matrix does not change. 186 * 187 * Example: \include SelfAdjointEigenSolver_compute_MatrixType.cpp 188 * Output: \verbinclude SelfAdjointEigenSolver_compute_MatrixType.out 189 * 190 * \sa SelfAdjointEigenSolver(const MatrixType&, int) 191 */ 192 SelfAdjointEigenSolver& compute(const MatrixType& matrix, int options = ComputeEigenvectors); 193 194 /** \brief Computes eigendecomposition of given matrix using a direct algorithm 195 * 196 * This is a variant of compute(const MatrixType&, int options) which 197 * directly solves the underlying polynomial equation. 198 * 199 * Currently only 3x3 matrices for which the sizes are known at compile time are supported (e.g., Matrix3d). 200 * 201 * This method is usually significantly faster than the QR algorithm 202 * but it might also be less accurate. It is also worth noting that 203 * for 3x3 matrices it involves trigonometric operations which are 204 * not necessarily available for all scalar types. 205 * 206 * \sa compute(const MatrixType&, int options) 207 */ 208 SelfAdjointEigenSolver& computeDirect(const MatrixType& matrix, int options = ComputeEigenvectors); 209 210 /** \brief Returns the eigenvectors of given matrix. 211 * 212 * \returns A const reference to the matrix whose columns are the eigenvectors. 213 * 214 * \pre The eigenvectors have been computed before. 215 * 216 * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 217 * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 218 * eigenvectors are normalized to have (Euclidean) norm equal to one. If 219 * this object was used to solve the eigenproblem for the selfadjoint 220 * matrix \f$ A \f$, then the matrix returned by this function is the 221 * matrix \f$ V \f$ in the eigendecomposition \f$ A = V D V^{-1} \f$. 222 * 223 * Example: \include SelfAdjointEigenSolver_eigenvectors.cpp 224 * Output: \verbinclude SelfAdjointEigenSolver_eigenvectors.out 225 * 226 * \sa eigenvalues() 227 */ 228 const MatrixType& eigenvectors() const 229 { 230 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 231 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 232 return m_eivec; 233 } 234 235 /** \brief Returns the eigenvalues of given matrix. 236 * 237 * \returns A const reference to the column vector containing the eigenvalues. 238 * 239 * \pre The eigenvalues have been computed before. 240 * 241 * The eigenvalues are repeated according to their algebraic multiplicity, 242 * so there are as many eigenvalues as rows in the matrix. The eigenvalues 243 * are sorted in increasing order. 244 * 245 * Example: \include SelfAdjointEigenSolver_eigenvalues.cpp 246 * Output: \verbinclude SelfAdjointEigenSolver_eigenvalues.out 247 * 248 * \sa eigenvectors(), MatrixBase::eigenvalues() 249 */ 250 const RealVectorType& eigenvalues() const 251 { 252 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 253 return m_eivalues; 254 } 255 256 /** \brief Computes the positive-definite square root of the matrix. 257 * 258 * \returns the positive-definite square root of the matrix 259 * 260 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 261 * have been computed before. 262 * 263 * The square root of a positive-definite matrix \f$ A \f$ is the 264 * positive-definite matrix whose square equals \f$ A \f$. This function 265 * uses the eigendecomposition \f$ A = V D V^{-1} \f$ to compute the 266 * square root as \f$ A^{1/2} = V D^{1/2} V^{-1} \f$. 267 * 268 * Example: \include SelfAdjointEigenSolver_operatorSqrt.cpp 269 * Output: \verbinclude SelfAdjointEigenSolver_operatorSqrt.out 270 * 271 * \sa operatorInverseSqrt(), 272 * \ref MatrixFunctions_Module "MatrixFunctions Module" 273 */ 274 MatrixType operatorSqrt() const 275 { 276 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 277 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 278 return m_eivec * m_eivalues.cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 279 } 280 281 /** \brief Computes the inverse square root of the matrix. 282 * 283 * \returns the inverse positive-definite square root of the matrix 284 * 285 * \pre The eigenvalues and eigenvectors of a positive-definite matrix 286 * have been computed before. 287 * 288 * This function uses the eigendecomposition \f$ A = V D V^{-1} \f$ to 289 * compute the inverse square root as \f$ V D^{-1/2} V^{-1} \f$. This is 290 * cheaper than first computing the square root with operatorSqrt() and 291 * then its inverse with MatrixBase::inverse(). 292 * 293 * Example: \include SelfAdjointEigenSolver_operatorInverseSqrt.cpp 294 * Output: \verbinclude SelfAdjointEigenSolver_operatorInverseSqrt.out 295 * 296 * \sa operatorSqrt(), MatrixBase::inverse(), 297 * \ref MatrixFunctions_Module "MatrixFunctions Module" 298 */ 299 MatrixType operatorInverseSqrt() const 300 { 301 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 302 eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 303 return m_eivec * m_eivalues.cwiseInverse().cwiseSqrt().asDiagonal() * m_eivec.adjoint(); 304 } 305 306 /** \brief Reports whether previous computation was successful. 307 * 308 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 309 */ 310 ComputationInfo info() const 311 { 312 eigen_assert(m_isInitialized && "SelfAdjointEigenSolver is not initialized."); 313 return m_info; 314 } 315 316 /** \brief Maximum number of iterations. 317 * 318 * The algorithm terminates if it does not converge within m_maxIterations * n iterations, where n 319 * denotes the size of the matrix. This value is currently set to 30 (copied from LAPACK). 320 */ 321 static const int m_maxIterations = 30; 322 323 #ifdef EIGEN2_SUPPORT 324 SelfAdjointEigenSolver(const MatrixType& matrix, bool computeEigenvectors) 325 : m_eivec(matrix.rows(), matrix.cols()), 326 m_eivalues(matrix.cols()), 327 m_subdiag(matrix.rows() > 1 ? matrix.rows() - 1 : 1), 328 m_isInitialized(false) 329 { 330 compute(matrix, computeEigenvectors); 331 } 332 333 SelfAdjointEigenSolver(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) 334 : m_eivec(matA.cols(), matA.cols()), 335 m_eivalues(matA.cols()), 336 m_subdiag(matA.cols() > 1 ? matA.cols() - 1 : 1), 337 m_isInitialized(false) 338 { 339 static_cast<GeneralizedSelfAdjointEigenSolver<MatrixType>*>(this)->compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); 340 } 341 342 void compute(const MatrixType& matrix, bool computeEigenvectors) 343 { 344 compute(matrix, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); 345 } 346 347 void compute(const MatrixType& matA, const MatrixType& matB, bool computeEigenvectors = true) 348 { 349 compute(matA, matB, computeEigenvectors ? ComputeEigenvectors : EigenvaluesOnly); 350 } 351 #endif // EIGEN2_SUPPORT 352 353 protected: 354 MatrixType m_eivec; 355 RealVectorType m_eivalues; 356 typename TridiagonalizationType::SubDiagonalType m_subdiag; 357 ComputationInfo m_info; 358 bool m_isInitialized; 359 bool m_eigenvectorsOk; 360 }; 361 362 /** \internal 363 * 364 * \eigenvalues_module \ingroup Eigenvalues_Module 365 * 366 * Performs a QR step on a tridiagonal symmetric matrix represented as a 367 * pair of two vectors \a diag and \a subdiag. 368 * 369 * \param matA the input selfadjoint matrix 370 * \param hCoeffs returned Householder coefficients 371 * 372 * For compilation efficiency reasons, this procedure does not use eigen expression 373 * for its arguments. 374 * 375 * Implemented from Golub's "Matrix Computations", algorithm 8.3.2: 376 * "implicit symmetric QR step with Wilkinson shift" 377 */ 378 namespace internal { 379 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> 380 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n); 381 } 382 383 template<typename MatrixType> 384 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> 385 ::compute(const MatrixType& matrix, int options) 386 { 387 using std::abs; 388 eigen_assert(matrix.cols() == matrix.rows()); 389 eigen_assert((options&~(EigVecMask|GenEigMask))==0 390 && (options&EigVecMask)!=EigVecMask 391 && "invalid option parameter"); 392 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 393 Index n = matrix.cols(); 394 m_eivalues.resize(n,1); 395 396 if(n==1) 397 { 398 m_eivalues.coeffRef(0,0) = numext::real(matrix.coeff(0,0)); 399 if(computeEigenvectors) 400 m_eivec.setOnes(n,n); 401 m_info = Success; 402 m_isInitialized = true; 403 m_eigenvectorsOk = computeEigenvectors; 404 return *this; 405 } 406 407 // declare some aliases 408 RealVectorType& diag = m_eivalues; 409 MatrixType& mat = m_eivec; 410 411 // map the matrix coefficients to [-1:1] to avoid over- and underflow. 412 mat = matrix.template triangularView<Lower>(); 413 RealScalar scale = mat.cwiseAbs().maxCoeff(); 414 if(scale==RealScalar(0)) scale = RealScalar(1); 415 mat.template triangularView<Lower>() /= scale; 416 m_subdiag.resize(n-1); 417 internal::tridiagonalization_inplace(mat, diag, m_subdiag, computeEigenvectors); 418 419 Index end = n-1; 420 Index start = 0; 421 Index iter = 0; // total number of iterations 422 423 while (end>0) 424 { 425 for (Index i = start; i<end; ++i) 426 if (internal::isMuchSmallerThan(abs(m_subdiag[i]),(abs(diag[i])+abs(diag[i+1])))) 427 m_subdiag[i] = 0; 428 429 // find the largest unreduced block 430 while (end>0 && m_subdiag[end-1]==0) 431 { 432 end--; 433 } 434 if (end<=0) 435 break; 436 437 // if we spent too many iterations, we give up 438 iter++; 439 if(iter > m_maxIterations * n) break; 440 441 start = end - 1; 442 while (start>0 && m_subdiag[start-1]!=0) 443 start--; 444 445 internal::tridiagonal_qr_step<MatrixType::Flags&RowMajorBit ? RowMajor : ColMajor>(diag.data(), m_subdiag.data(), start, end, computeEigenvectors ? m_eivec.data() : (Scalar*)0, n); 446 } 447 448 if (iter <= m_maxIterations * n) 449 m_info = Success; 450 else 451 m_info = NoConvergence; 452 453 // Sort eigenvalues and corresponding vectors. 454 // TODO make the sort optional ? 455 // TODO use a better sort algorithm !! 456 if (m_info == Success) 457 { 458 for (Index i = 0; i < n-1; ++i) 459 { 460 Index k; 461 m_eivalues.segment(i,n-i).minCoeff(&k); 462 if (k > 0) 463 { 464 std::swap(m_eivalues[i], m_eivalues[k+i]); 465 if(computeEigenvectors) 466 m_eivec.col(i).swap(m_eivec.col(k+i)); 467 } 468 } 469 } 470 471 // scale back the eigen values 472 m_eivalues *= scale; 473 474 m_isInitialized = true; 475 m_eigenvectorsOk = computeEigenvectors; 476 return *this; 477 } 478 479 480 namespace internal { 481 482 template<typename SolverType,int Size,bool IsComplex> struct direct_selfadjoint_eigenvalues 483 { 484 static inline void run(SolverType& eig, const typename SolverType::MatrixType& A, int options) 485 { eig.compute(A,options); } 486 }; 487 488 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,3,false> 489 { 490 typedef typename SolverType::MatrixType MatrixType; 491 typedef typename SolverType::RealVectorType VectorType; 492 typedef typename SolverType::Scalar Scalar; 493 494 static inline void computeRoots(const MatrixType& m, VectorType& roots) 495 { 496 using std::sqrt; 497 using std::atan2; 498 using std::cos; 499 using std::sin; 500 const Scalar s_inv3 = Scalar(1.0)/Scalar(3.0); 501 const Scalar s_sqrt3 = sqrt(Scalar(3.0)); 502 503 // The characteristic equation is x^3 - c2*x^2 + c1*x - c0 = 0. The 504 // eigenvalues are the roots to this equation, all guaranteed to be 505 // real-valued, because the matrix is symmetric. 506 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); 507 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); 508 Scalar c2 = m(0,0) + m(1,1) + m(2,2); 509 510 // Construct the parameters used in classifying the roots of the equation 511 // and in solving the equation for the roots in closed form. 512 Scalar c2_over_3 = c2*s_inv3; 513 Scalar a_over_3 = (c1 - c2*c2_over_3)*s_inv3; 514 if (a_over_3 > Scalar(0)) 515 a_over_3 = Scalar(0); 516 517 Scalar half_b = Scalar(0.5)*(c0 + c2_over_3*(Scalar(2)*c2_over_3*c2_over_3 - c1)); 518 519 Scalar q = half_b*half_b + a_over_3*a_over_3*a_over_3; 520 if (q > Scalar(0)) 521 q = Scalar(0); 522 523 // Compute the eigenvalues by solving for the roots of the polynomial. 524 Scalar rho = sqrt(-a_over_3); 525 Scalar theta = atan2(sqrt(-q),half_b)*s_inv3; 526 Scalar cos_theta = cos(theta); 527 Scalar sin_theta = sin(theta); 528 roots(0) = c2_over_3 + Scalar(2)*rho*cos_theta; 529 roots(1) = c2_over_3 - rho*(cos_theta + s_sqrt3*sin_theta); 530 roots(2) = c2_over_3 - rho*(cos_theta - s_sqrt3*sin_theta); 531 532 // Sort in increasing order. 533 if (roots(0) >= roots(1)) 534 std::swap(roots(0),roots(1)); 535 if (roots(1) >= roots(2)) 536 { 537 std::swap(roots(1),roots(2)); 538 if (roots(0) >= roots(1)) 539 std::swap(roots(0),roots(1)); 540 } 541 } 542 543 static inline void run(SolverType& solver, const MatrixType& mat, int options) 544 { 545 using std::sqrt; 546 eigen_assert(mat.cols() == 3 && mat.cols() == mat.rows()); 547 eigen_assert((options&~(EigVecMask|GenEigMask))==0 548 && (options&EigVecMask)!=EigVecMask 549 && "invalid option parameter"); 550 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 551 552 MatrixType& eivecs = solver.m_eivec; 553 VectorType& eivals = solver.m_eivalues; 554 555 // map the matrix coefficients to [-1:1] to avoid over- and underflow. 556 Scalar scale = mat.cwiseAbs().maxCoeff(); 557 MatrixType scaledMat = mat / scale; 558 559 // compute the eigenvalues 560 computeRoots(scaledMat,eivals); 561 562 // compute the eigen vectors 563 if(computeEigenvectors) 564 { 565 Scalar safeNorm2 = Eigen::NumTraits<Scalar>::epsilon(); 566 safeNorm2 *= safeNorm2; 567 if((eivals(2)-eivals(0))<=Eigen::NumTraits<Scalar>::epsilon()) 568 { 569 eivecs.setIdentity(); 570 } 571 else 572 { 573 scaledMat = scaledMat.template selfadjointView<Lower>(); 574 MatrixType tmp; 575 tmp = scaledMat; 576 577 Scalar d0 = eivals(2) - eivals(1); 578 Scalar d1 = eivals(1) - eivals(0); 579 int k = d0 > d1 ? 2 : 0; 580 d0 = d0 > d1 ? d1 : d0; 581 582 tmp.diagonal().array () -= eivals(k); 583 VectorType cross; 584 Scalar n; 585 n = (cross = tmp.row(0).cross(tmp.row(1))).squaredNorm(); 586 587 if(n>safeNorm2) 588 eivecs.col(k) = cross / sqrt(n); 589 else 590 { 591 n = (cross = tmp.row(0).cross(tmp.row(2))).squaredNorm(); 592 593 if(n>safeNorm2) 594 eivecs.col(k) = cross / sqrt(n); 595 else 596 { 597 n = (cross = tmp.row(1).cross(tmp.row(2))).squaredNorm(); 598 599 if(n>safeNorm2) 600 eivecs.col(k) = cross / sqrt(n); 601 else 602 { 603 // the input matrix and/or the eigenvaues probably contains some inf/NaN, 604 // => exit 605 // scale back to the original size. 606 eivals *= scale; 607 608 solver.m_info = NumericalIssue; 609 solver.m_isInitialized = true; 610 solver.m_eigenvectorsOk = computeEigenvectors; 611 return; 612 } 613 } 614 } 615 616 tmp = scaledMat; 617 tmp.diagonal().array() -= eivals(1); 618 619 if(d0<=Eigen::NumTraits<Scalar>::epsilon()) 620 eivecs.col(1) = eivecs.col(k).unitOrthogonal(); 621 else 622 { 623 n = (cross = eivecs.col(k).cross(tmp.row(0).normalized())).squaredNorm(); 624 if(n>safeNorm2) 625 eivecs.col(1) = cross / sqrt(n); 626 else 627 { 628 n = (cross = eivecs.col(k).cross(tmp.row(1))).squaredNorm(); 629 if(n>safeNorm2) 630 eivecs.col(1) = cross / sqrt(n); 631 else 632 { 633 n = (cross = eivecs.col(k).cross(tmp.row(2))).squaredNorm(); 634 if(n>safeNorm2) 635 eivecs.col(1) = cross / sqrt(n); 636 else 637 { 638 // we should never reach this point, 639 // if so the last two eigenvalues are likely to ve very closed to each other 640 eivecs.col(1) = eivecs.col(k).unitOrthogonal(); 641 } 642 } 643 } 644 645 // make sure that eivecs[1] is orthogonal to eivecs[2] 646 Scalar d = eivecs.col(1).dot(eivecs.col(k)); 647 eivecs.col(1) = (eivecs.col(1) - d * eivecs.col(k)).normalized(); 648 } 649 650 eivecs.col(k==2 ? 0 : 2) = eivecs.col(k).cross(eivecs.col(1)).normalized(); 651 } 652 } 653 // Rescale back to the original size. 654 eivals *= scale; 655 656 solver.m_info = Success; 657 solver.m_isInitialized = true; 658 solver.m_eigenvectorsOk = computeEigenvectors; 659 } 660 }; 661 662 // 2x2 direct eigenvalues decomposition, code from Hauke Heibel 663 template<typename SolverType> struct direct_selfadjoint_eigenvalues<SolverType,2,false> 664 { 665 typedef typename SolverType::MatrixType MatrixType; 666 typedef typename SolverType::RealVectorType VectorType; 667 typedef typename SolverType::Scalar Scalar; 668 669 static inline void computeRoots(const MatrixType& m, VectorType& roots) 670 { 671 using std::sqrt; 672 const Scalar t0 = Scalar(0.5) * sqrt( numext::abs2(m(0,0)-m(1,1)) + Scalar(4)*m(1,0)*m(1,0)); 673 const Scalar t1 = Scalar(0.5) * (m(0,0) + m(1,1)); 674 roots(0) = t1 - t0; 675 roots(1) = t1 + t0; 676 } 677 678 static inline void run(SolverType& solver, const MatrixType& mat, int options) 679 { 680 using std::sqrt; 681 eigen_assert(mat.cols() == 2 && mat.cols() == mat.rows()); 682 eigen_assert((options&~(EigVecMask|GenEigMask))==0 683 && (options&EigVecMask)!=EigVecMask 684 && "invalid option parameter"); 685 bool computeEigenvectors = (options&ComputeEigenvectors)==ComputeEigenvectors; 686 687 MatrixType& eivecs = solver.m_eivec; 688 VectorType& eivals = solver.m_eivalues; 689 690 // map the matrix coefficients to [-1:1] to avoid over- and underflow. 691 Scalar scale = mat.cwiseAbs().maxCoeff(); 692 scale = (std::max)(scale,Scalar(1)); 693 MatrixType scaledMat = mat / scale; 694 695 // Compute the eigenvalues 696 computeRoots(scaledMat,eivals); 697 698 // compute the eigen vectors 699 if(computeEigenvectors) 700 { 701 scaledMat.diagonal().array () -= eivals(1); 702 Scalar a2 = numext::abs2(scaledMat(0,0)); 703 Scalar c2 = numext::abs2(scaledMat(1,1)); 704 Scalar b2 = numext::abs2(scaledMat(1,0)); 705 if(a2>c2) 706 { 707 eivecs.col(1) << -scaledMat(1,0), scaledMat(0,0); 708 eivecs.col(1) /= sqrt(a2+b2); 709 } 710 else 711 { 712 eivecs.col(1) << -scaledMat(1,1), scaledMat(1,0); 713 eivecs.col(1) /= sqrt(c2+b2); 714 } 715 716 eivecs.col(0) << eivecs.col(1).unitOrthogonal(); 717 } 718 719 // Rescale back to the original size. 720 eivals *= scale; 721 722 solver.m_info = Success; 723 solver.m_isInitialized = true; 724 solver.m_eigenvectorsOk = computeEigenvectors; 725 } 726 }; 727 728 } 729 730 template<typename MatrixType> 731 SelfAdjointEigenSolver<MatrixType>& SelfAdjointEigenSolver<MatrixType> 732 ::computeDirect(const MatrixType& matrix, int options) 733 { 734 internal::direct_selfadjoint_eigenvalues<SelfAdjointEigenSolver,Size,NumTraits<Scalar>::IsComplex>::run(*this,matrix,options); 735 return *this; 736 } 737 738 namespace internal { 739 template<int StorageOrder,typename RealScalar, typename Scalar, typename Index> 740 static void tridiagonal_qr_step(RealScalar* diag, RealScalar* subdiag, Index start, Index end, Scalar* matrixQ, Index n) 741 { 742 using std::abs; 743 RealScalar td = (diag[end-1] - diag[end])*RealScalar(0.5); 744 RealScalar e = subdiag[end-1]; 745 // Note that thanks to scaling, e^2 or td^2 cannot overflow, however they can still 746 // underflow thus leading to inf/NaN values when using the following commented code: 747 // RealScalar e2 = numext::abs2(subdiag[end-1]); 748 // RealScalar mu = diag[end] - e2 / (td + (td>0 ? 1 : -1) * sqrt(td*td + e2)); 749 // This explain the following, somewhat more complicated, version: 750 RealScalar mu = diag[end]; 751 if(td==0) 752 mu -= abs(e); 753 else 754 { 755 RealScalar e2 = numext::abs2(subdiag[end-1]); 756 RealScalar h = numext::hypot(td,e); 757 if(e2==0) mu -= (e / (td + (td>0 ? 1 : -1))) * (e / h); 758 else mu -= e2 / (td + (td>0 ? h : -h)); 759 } 760 761 RealScalar x = diag[start] - mu; 762 RealScalar z = subdiag[start]; 763 for (Index k = start; k < end; ++k) 764 { 765 JacobiRotation<RealScalar> rot; 766 rot.makeGivens(x, z); 767 768 // do T = G' T G 769 RealScalar sdk = rot.s() * diag[k] + rot.c() * subdiag[k]; 770 RealScalar dkp1 = rot.s() * subdiag[k] + rot.c() * diag[k+1]; 771 772 diag[k] = rot.c() * (rot.c() * diag[k] - rot.s() * subdiag[k]) - rot.s() * (rot.c() * subdiag[k] - rot.s() * diag[k+1]); 773 diag[k+1] = rot.s() * sdk + rot.c() * dkp1; 774 subdiag[k] = rot.c() * sdk - rot.s() * dkp1; 775 776 777 if (k > start) 778 subdiag[k - 1] = rot.c() * subdiag[k-1] - rot.s() * z; 779 780 x = subdiag[k]; 781 782 if (k < end - 1) 783 { 784 z = -rot.s() * subdiag[k+1]; 785 subdiag[k + 1] = rot.c() * subdiag[k+1]; 786 } 787 788 // apply the givens rotation to the unit matrix Q = Q * G 789 if (matrixQ) 790 { 791 // FIXME if StorageOrder == RowMajor this operation is not very efficient 792 Map<Matrix<Scalar,Dynamic,Dynamic,StorageOrder> > q(matrixQ,n,n); 793 q.applyOnTheRight(k,k+1,rot); 794 } 795 } 796 } 797 798 } // end namespace internal 799 800 } // end namespace Eigen 801 802 #endif // EIGEN_SELFADJOINTEIGENSOLVER_H 803