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