EigenSolver.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
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> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 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 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Sets the maximum number of iterations allowed. */ 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EigenSolver& setMaxIterations(Index maxIters) 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realSchur.setMaxIterations(maxIters); 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the maximum number of iterations. */ 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index getMaxIterations() 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_realSchur.getMaxIterations(); 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void doComputeEigenvectors(); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_eivec; 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvalueType m_eivalues; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_eigenvectorsOk; 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> m_realSchur; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matT; 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColumnVectorType m_tmp; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivalues.rows(); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matD = MatrixType::Zero(n,n); 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)))) 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i)); 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)), 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)); 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matD; 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivec.cols(); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorsType matV(n,n); 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<n; ++j) 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j))) || j+1==n) 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a real eigen value 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a pair of complex eigen values 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j+1).normalize(); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++j; 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matV; 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezEigenSolver<MatrixType>& 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(matrix.cols() == matrix.rows()); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reduce to real Schur form. 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_realSchur.compute(matrix, computeEigenvectors); 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_realSchur.info() == Success) 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT = m_realSchur.matrixT(); 377c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 378c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec = m_realSchur.matrixU(); 379c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvalues from matT 381c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.resize(matrix.cols()); 382c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = 0; 383c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (i < matrix.cols()) 384c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 385c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 386c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = m_matT.coeff(i, i); 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 389c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 390c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 391c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 392c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath i += 2; 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvectors. 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doComputeEigenvectors(); 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk = computeEigenvectors; 407c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 408c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 409c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 410c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 411c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Complex scalar division. 412c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename Scalar> 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandezstd::complex<Scalar> cdiv(const Scalar& xr, const Scalar& xi, const Scalar& yr, const Scalar& yi) 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r,d; 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(yr) > abs(yi)) 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 419c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = yi/yr; 420c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d = yr + r*yi; 421c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return std::complex<Scalar>((xr + r*xi)/d, (xi - r*xr)/d); 422c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 423c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 424c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 425c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath r = yr/yi; 426c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath d = yi + r*yr; 427c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return std::complex<Scalar>((r*xr + xi)/d, (r*xi - xr)/d); 428c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 429c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 430c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 431c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid EigenSolver<MatrixType>::doComputeEigenvectors() 434c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 436c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_eivec.cols(); 437c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar eps = NumTraits<Scalar>::epsilon(); 438c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 439c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // inefficient! this is already computed in RealSchur 440c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar norm(0); 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < size; ++j) 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Backsubstitute to find vectors of upper triangular form 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (norm == 0.0) 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index n = size-1; n >= 0; n--) 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = m_eivalues.coeff(n).real(); 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar q = m_eivalues.coeff(n).imag(); 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Scalar vector 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (q == Scalar(0)) 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 460c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastr(0), lastw(0); 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n; 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n) = 1.0; 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-1; i >= 0; i--) 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() < 0.0) 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 472c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastr = r; 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() == 0.0) 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (w != 0.0) 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / w; 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / (eps * norm); 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else // Solve real equations 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 488c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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(); 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = (x * lastr - lastw * r) / denom; 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = t; 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(x) > abs(lastw)) 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-r - w * t) / x; 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 494c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t = abs(m_matT.coeff(i,n)); 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.col(n).tail(size-i) /= t; 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 502c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 504c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (q < Scalar(0) && n > 0) // Complex vector 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastra(0), lastsa(0), lastw(0); 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n-1; 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Last vector component imaginary so matrix is triangular 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n))) 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 516c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::complex<Scalar> cc = cdiv<Scalar>(0.0,-m_matT.coeff(n-1,n),m_matT.coeff(n-1,n-1)-p,q); 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(n-1,n-1) = numext::real(cc); 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(n-1,n) = numext::imag(cc); 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n-1) = 0.0; 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n,n) = 1.0; 523c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-2; i >= 0; i--) 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() < 0.0) 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastra = ra; 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastsa = sa; 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 535c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() == RealScalar(0)) 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::complex<Scalar> cc = cdiv(-ra,-sa,w,q); 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n-1) = numext::real(cc); 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n) = numext::imag(cc); 543c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 544c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 546c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Solve complex equations 547c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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; 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((vr == 0.0) && (vi == 0.0)) 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez std::complex<Scalar> cc = cdiv(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra,vr,vi); 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n-1) = numext::real(cc); 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n) = numext::imag(cc); 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(x) > (abs(lastw) + abs(q))) 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath cc = cdiv(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n),lastw,q); 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i+1,n-1) = numext::real(cc); 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i+1,n) = numext::imag(cc); 567c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath using std::max; 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t = (max)(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.block(i, n-1, size-i, 2) /= t; 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 576c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 577c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 579c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // We handled a pair of complex conjugate eigenvalues, so need to skip them both 580c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n--; 581c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 582c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(0 && "Internal bug in EigenSolver"); // this should not happen 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Back transformation to get eigenvectors of original matrix 589c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = size-1; j >= 0; j--) 590c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 591c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1); 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(j) = m_tmp; 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 596c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_EIGENSOLVER_H 599