EigenSolver.h revision c981c48f5bc9aefeffc0bcb0cc3934c2fae179dd
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 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_EIGENSOLVER_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_EIGENSOLVER_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "./RealSchur.h" 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class EigenSolver 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes eigenvalues and eigenvectors of general matrices 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are computing the 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigendecomposition; this is expected to be an instantiation of the Matrix 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * class template. Currently, only real matrices are supported. 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v \f$. If 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * have \f$ A = V D V^{-1} \f$. This is called the eigendecomposition. 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues and eigenvectors of a matrix may be complex, even when the 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix is real. However, we can choose real matrices \f$ V \f$ and \f$ D 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ satisfying \f$ A V = V D \f$, just like the eigendecomposition, if the 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix \f$ D \f$ is not required to be diagonal, but if it is allowed to 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * have blocks of the form 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f[ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f] 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (where \f$ u \f$ and \f$ v \f$ are real numbers) on the diagonal. These 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * blocks correspond to complex eigenvalue pairs \f$ u \pm iv \f$. We call 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this variant of the eigendecomposition the pseudo-eigendecomposition. 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Call the function compute() to compute the eigenvalues and eigenvectors of 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * a given matrix. Alternatively, you can use the 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&, bool) constructor which computes the 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues and eigenvectors at construction time. Once the eigenvalue and 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors are computed, they can be retrieved with the eigenvalues() and 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors() functions. The pseudoEigenvalueMatrix() and 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * pseudoEigenvectors() methods allow the construction of the 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * pseudo-eigendecomposition. 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The documentation for EigenSolver(const MatrixType&, bool) contains an 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * example of the typical use of this class. 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \note The implementation is adapted from 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * <a href="http://math.nist.gov/javanumerics/jama/">JAMA</a> (public domain). 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Their code is based on EISPACK. 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class EigenSolver 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Synonym for the template parameter \p _MatrixType. */ 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Scalar type for matrices of type #MatrixType. */ 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Complex scalar type for #MatrixType. 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is \c std::complex<Scalar> if #Scalar is real (e.g., 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c float or \c double) and just \c Scalar if #Scalar is 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * complex. 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::complex<RealScalar> ComplexScalar; 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a column vector with entries of type #ComplexScalar. 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The length of the vector is the size of #MatrixType. 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a square matrix with entries of type #ComplexScalar. 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The size is the same as the size of #MatrixType. 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor. 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via EigenSolver::compute(const MatrixType&, bool). 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example. 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver() : m_eivec(), m_eivalues(), m_isInitialized(false), m_realSchur(), m_matT(), m_tmp() {} 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor with memory preallocation 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Like the default constructor but with preallocation of the internal data 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * according to the specified problem \a size. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa EigenSolver() 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver(Index size) 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(size, size), 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(size), 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_realSchur(size), 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT(size, size), 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp(size) 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes eigendecomposition of given matrix. 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute() to compute the eigenvalues 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and eigenvectors. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_EigenSolver_MatrixType.cpp 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_EigenSolver_MatrixType.out 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(matrix.rows(), matrix.cols()), 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(matrix.cols()), 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_realSchur(matrix.cols()), 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT(matrix.rows(), matrix.cols()), 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp(matrix.cols()) 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix, computeEigenvectors); 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvectors of given matrix. 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns %Matrix whose columns are the (possibly complex) eigenvectors. 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before, and 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors are normalized to have (Euclidean) norm equal to one. The 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix returned by this function is the matrix \f$ V \f$ in the 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists. 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_eigenvectors.cpp 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_eigenvectors.out 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa eigenvalues(), pseudoEigenvectors() 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorsType eigenvectors() const; 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the pseudo-eigenvectors of given matrix. 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Const reference to matrix whose columns are the pseudo-eigenvectors. 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before, and 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The real matrix \f$ V \f$ returned by this function and the 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix() 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * satisfy \f$ AV = VD \f$. 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_pseudoEigenvectors.cpp 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_pseudoEigenvectors.out 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa pseudoEigenvalueMatrix(), eigenvectors() 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& pseudoEigenvectors() const 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivec; 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition. 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A block-diagonal matrix. 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before. 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix \f$ D \f$ returned by this function is real and 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * blocks of the form 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$. 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * These blocks are not sorted in any particular order. 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * pseudoEigenvectors() satisfy \f$ AV = VD \f$. 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa pseudoEigenvectors() for an example, eigenvalues() 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType pseudoEigenvalueMatrix() const; 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvalues of given matrix. 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the column vector containing the eigenvalues. 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before. 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues are repeated according to their algebraic multiplicity, 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * so there are as many eigenvalues as rows in the matrix. The eigenvalues 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * are not sorted in any particular order. 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_eigenvalues.cpp 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_eigenvalues.out 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa eigenvectors(), pseudoEigenvalueMatrix(), 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::eigenvalues() 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvalueType& eigenvalues() const 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivalues; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes eigendecomposition of given matrix. 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function computes the eigenvalues of the real matrix \p matrix. 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve them. If 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors is true, then the eigenvectors are also computed 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and can be retrieved by calling eigenvectors(). 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix is first reduced to real Schur form using the RealSchur 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * class. The Schur decomposition is then used to compute the eigenvalues 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and eigenvectors. 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The cost of the computation is dominated by the cost of the 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Schur decomposition, which is very approximately \f$ 25n^3 \f$ 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false. 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method reuses of the allocated data in the EigenSolver object. 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_compute.cpp 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_compute.out 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_realSchur.info(); 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void doComputeEigenvectors(); 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_eivec; 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvalueType m_eivalues; 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_eigenvectorsOk; 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> m_realSchur; 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matT; 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColumnVectorType m_tmp; 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivalues.rows(); 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matD = MatrixType::Zero(n,n); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i)))) 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matD.coeffRef(i,i) = internal::real(m_eivalues.coeff(i)); 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matD.template block<2,2>(i,i) << internal::real(m_eivalues.coeff(i)), internal::imag(m_eivalues.coeff(i)), 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath -internal::imag(m_eivalues.coeff(i)), internal::real(m_eivalues.coeff(i)); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matD; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivec.cols(); 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorsType matV(n,n); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<n; ++j) 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::isMuchSmallerThan(internal::imag(m_eivalues.coeff(j)), internal::real(m_eivalues.coeff(j))) || j+1==n) 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a real eigen value 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a pair of complex eigen values 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j+1).normalize(); 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++j; 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matV; 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathEigenSolver<MatrixType>& EigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(matrix.cols() == matrix.rows()); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reduce to real Schur form. 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_realSchur.compute(matrix, computeEigenvectors); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_realSchur.info() == Success) 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT = m_realSchur.matrixT(); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec = m_realSchur.matrixU(); 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvalues from matT 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.resize(matrix.cols()); 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = 0; 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (i < matrix.cols()) 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = m_matT.coeff(i, i); 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar z = internal::sqrt(internal::abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath i += 2; 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvectors. 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doComputeEigenvectors(); 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk = computeEigenvectors; 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 393c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Complex scalar division. 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstd::complex<Scalar> cdiv(Scalar xr, Scalar xi, Scalar yr, Scalar yi) 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r,d; 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(yr) > internal::abs(yi)) 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = yi/yr; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d = yr + r*yi; 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d); 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = yr/yi; 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d = yi + r*yr; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d); 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 413c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid EigenSolver<MatrixType>::doComputeEigenvectors() 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_eivec.cols(); 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar eps = NumTraits<Scalar>::epsilon(); 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // inefficient! this is already computed in RealSchur 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar norm(0); 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < size; ++j) 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Backsubstitute to find vectors of upper triangular form 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (norm == 0.0) 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index n = size-1; n >= 0; n--) 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 435c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = m_eivalues.coeff(n).real(); 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar q = m_eivalues.coeff(n).imag(); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Scalar vector 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (q == Scalar(0)) 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastr(0), lastw(0); 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n; 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n) = 1.0; 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-1; i >= 0; i--) 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() < 0.0) 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastr = r; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() == 0.0) 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (w != 0.0) 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / w; 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / (eps * norm); 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else // Solve real equations 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar denom = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag(); 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = (x * lastr - lastw * r) / denom; 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = t; 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(x) > internal::abs(lastw)) 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-r - w * t) / x; 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = internal::abs(m_matT.coeff(i,n)); 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.col(n).tail(size-i) /= t; 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (q < Scalar(0) && n > 0) // Complex vector 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastra(0), lastsa(0), lastw(0); 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n-1; 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Last vector component imaginary so matrix is triangular 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(m_matT.coeff(n,n-1)) > internal::abs(m_matT.coeff(n-1,n))) 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n-1) = internal::real(cc); 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n) = internal::imag(cc); 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n-1) = 0.0; 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n) = 1.0; 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-2; i >= 0; i--) 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() < 0.0) 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastra = ra; 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastsa = sa; 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() == RealScalar(0)) 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n-1) = internal::real(cc); 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = internal::imag(cc); 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Solve complex equations 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar vr = (m_eivalues.coeff(i).real() - p) * (m_eivalues.coeff(i).real() - p) + m_eivalues.coeff(i).imag() * m_eivalues.coeff(i).imag() - q * q; 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((vr == 0.0) && (vi == 0.0)) 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath vr = eps * norm * (internal::abs(w) + internal::abs(q) + internal::abs(x) + internal::abs(y) + internal::abs(lastw)); 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n-1) = internal::real(cc); 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = internal::imag(cc); 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (internal::abs(x) > (internal::abs(lastw) + internal::abs(q))) 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; 542c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n-1) = internal::real(cc); 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = internal::imag(cc); 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::max; 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = (max)(internal::abs(m_matT.coeff(i,n-1)),internal::abs(m_matT.coeff(i,n))); 554c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.block(i, n-1, size-i, 2) /= t; 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // We handled a pair of complex conjugate eigenvalues, so need to skip them both 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n--; 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 565c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen 566c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Back transformation to get eigenvectors of original matrix 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = size-1; j >= 0; j--) 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(j) = m_tmp; 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_EIGENSOLVER_H 580