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; 822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 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 */ 1132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 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 */ 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit 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 */ 1462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit EigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(matrix.rows(), matrix.cols()), 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(matrix.cols()), 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_realSchur(matrix.cols()), 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT(matrix.rows(), matrix.cols()), 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp(matrix.cols()) 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang compute(matrix.derived(), computeEigenvectors); 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvectors of given matrix. 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns %Matrix whose columns are the (possibly complex) eigenvectors. 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before, and 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Column \f$ k \f$ of the returned matrix is an eigenvector corresponding 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to eigenvalue number \f$ k \f$ as returned by eigenvalues(). The 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvectors are normalized to have (Euclidean) norm equal to one. The 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix returned by this function is the matrix \f$ V \f$ in the 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigendecomposition \f$ A = V D V^{-1} \f$, if it exists. 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_eigenvectors.cpp 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_eigenvectors.out 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa eigenvalues(), pseudoEigenvectors() 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorsType eigenvectors() const; 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the pseudo-eigenvectors of given matrix. 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Const reference to matrix whose columns are the pseudo-eigenvectors. 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before, and 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The real matrix \f$ V \f$ returned by this function and the 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * block-diagonal matrix \f$ D \f$ returned by pseudoEigenvalueMatrix() 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * satisfy \f$ AV = VD \f$. 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_pseudoEigenvectors.cpp 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_pseudoEigenvectors.out 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa pseudoEigenvalueMatrix(), eigenvectors() 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& pseudoEigenvectors() const 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivec; 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the block-diagonal matrix in the pseudo-eigendecomposition. 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A block-diagonal matrix. 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before. 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix \f$ D \f$ returned by this function is real and 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * block-diagonal. The blocks on the diagonal are either 1-by-1 or 2-by-2 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * blocks of the form 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \begin{bmatrix} u & v \\ -v & u \end{bmatrix} \f$. 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * These blocks are not sorted in any particular order. 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix \f$ D \f$ and the matrix \f$ V \f$ returned by 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * pseudoEigenvectors() satisfy \f$ AV = VD \f$. 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa pseudoEigenvectors() for an example, eigenvalues() 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType pseudoEigenvalueMatrix() const; 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvalues of given matrix. 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the column vector containing the eigenvalues. 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * EigenSolver(const MatrixType&,bool) or the member function 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute(const MatrixType&, bool) has been called before. 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues are repeated according to their algebraic multiplicity, 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * so there are as many eigenvalues as rows in the matrix. The eigenvalues 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * are not sorted in any particular order. 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_eigenvalues.cpp 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_eigenvalues.out 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa eigenvectors(), pseudoEigenvalueMatrix(), 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * MatrixBase::eigenvalues() 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvalueType& eigenvalues() const 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivalues; 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes eigendecomposition of given matrix. 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function computes the eigenvalues of the real matrix \p matrix. 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve them. If 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors is true, then the eigenvectors are also computed 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and can be retrieved by calling eigenvectors(). 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix is first reduced to real Schur form using the RealSchur 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * class. The Schur decomposition is then used to compute the eigenvalues 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and eigenvectors. 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The cost of the computation is dominated by the cost of the 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Schur decomposition, which is very approximately \f$ 25n^3 \f$ 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (where \f$ n \f$ is the size of the matrix) if \p computeEigenvectors 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is true, and \f$ 10n^3 \f$ if \p computeEigenvectors is false. 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method reuses of the allocated data in the EigenSolver object. 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include EigenSolver_compute.cpp 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude EigenSolver_compute.out 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 2782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \returns NumericalIssue if the input contains INF or NaN values or overflow occured. Returns Success otherwise. */ 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 2842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_info; 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Sets the maximum number of iterations allowed. */ 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EigenSolver& setMaxIterations(Index maxIters) 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realSchur.setMaxIterations(maxIters); 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the maximum number of iterations. */ 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index getMaxIterations() 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_realSchur.getMaxIterations(); 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void doComputeEigenvectors(); 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 304a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 305a829215e078ace896f52702caa0c27608f40e3b0Miao Wang static void check_template_parameters() 306a829215e078ace896f52702caa0c27608f40e3b0Miao Wang { 307a829215e078ace896f52702caa0c27608f40e3b0Miao Wang EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 308a829215e078ace896f52702caa0c27608f40e3b0Miao Wang EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 309a829215e078ace896f52702caa0c27608f40e3b0Miao Wang } 310a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_eivec; 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvalueType m_eivalues; 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_eigenvectorsOk; 3152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComputationInfo m_info; 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealSchur<MatrixType> m_realSchur; 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matT; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColumnVectorType m_tmp; 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathMatrixType EigenSolver<MatrixType>::pseudoEigenvalueMatrix() const 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivalues.rows(); 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType matD = MatrixType::Zero(n,n); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)), precision)) 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez matD.coeffRef(i,i) = numext::real(m_eivalues.coeff(i)); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez matD.template block<2,2>(i,i) << numext::real(m_eivalues.coeff(i)), numext::imag(m_eivalues.coeff(i)), 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez -numext::imag(m_eivalues.coeff(i)), numext::real(m_eivalues.coeff(i)); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matD; 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtypename EigenSolver<MatrixType>::EigenvectorsType EigenSolver<MatrixType>::eigenvectors() const 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "EigenSolver is not initialized."); 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 3492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar precision = RealScalar(2)*NumTraits<RealScalar>::epsilon(); 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = m_eivec.cols(); 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorsType matV(n,n); 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j=0; j<n; ++j) 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 3542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (internal::isMuchSmallerThan(numext::imag(m_eivalues.coeff(j)), numext::real(m_eivalues.coeff(j)), precision) || j+1==n) 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a real eigen value 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j) = m_eivec.col(j).template cast<ComplexScalar>(); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // we have a pair of complex eigen values 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; ++i) 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j) = ComplexScalar(m_eivec.coeff(i,j), m_eivec.coeff(i,j+1)); 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.coeffRef(i,j+1) = ComplexScalar(m_eivec.coeff(i,j), -m_eivec.coeff(i,j+1)); 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j).normalize(); 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matV.col(j+1).normalize(); 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++j; 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return matV; 374c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 375c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 376c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename InputType> 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezEigenSolver<MatrixType>& 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) 380c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 381a829215e078ace896f52702caa0c27608f40e3b0Miao Wang check_template_parameters(); 382a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 3852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang using numext::isfinite; 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(matrix.cols() == matrix.rows()); 387c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 388c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Reduce to real Schur form. 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_realSchur.compute(matrix.derived(), computeEigenvectors); 3902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = m_realSchur.info(); 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_info == Success) 394c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 395c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT = m_realSchur.matrixT(); 396c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 397c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec = m_realSchur.matrixU(); 398c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 399c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvalues from matT 400c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.resize(matrix.cols()); 401c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index i = 0; 402c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath while (i < matrix.cols()) 403c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 404c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (i == matrix.cols() - 1 || m_matT.coeff(i+1, i) == Scalar(0)) 405c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 406c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = m_matT.coeff(i, i); 4072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!(isfinite)(m_eivalues.coeffRef(i))) 4082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = true; 4102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eigenvectorsOk = false; 4112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = NumericalIssue; 4122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *this; 4132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 414c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ++i; 415c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 416c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 417c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 418c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = Scalar(0.5) * (m_matT.coeff(i, i) - m_matT.coeff(i+1, i+1)); 4192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar z; 4202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Compute z = sqrt(abs(p * p + m_matT.coeff(i+1, i) * m_matT.coeff(i, i+1))); 4212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // without overflow 4222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar t0 = m_matT.coeff(i+1, i); 4242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar t1 = m_matT.coeff(i, i+1); 4252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar maxval = numext::maxi<Scalar>(abs(p),numext::maxi<Scalar>(abs(t0),abs(t1))); 4262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang t0 /= maxval; 4272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang t1 /= maxval; 4282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar p0 = p/maxval; 4292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang z = maxval * sqrt(abs(p0 * p0 + t0 * t1)); 4302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 432c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, z); 433c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.coeffRef(i+1) = ComplexScalar(m_matT.coeff(i+1, i+1) + p, -z); 4342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(!((isfinite)(m_eivalues.coeffRef(i)) && (isfinite)(m_eivalues.coeffRef(i+1)))) 4352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 4362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_isInitialized = true; 4372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eigenvectorsOk = false; 4382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_info = NumericalIssue; 4392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return *this; 4402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 441c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath i += 2; 442c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 443c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 444c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 445c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute eigenvectors. 446c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (computeEigenvectors) 447c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doComputeEigenvectors(); 448c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 449c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 450c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 451c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk = computeEigenvectors; 452c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 453c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 454c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 455c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 456c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 457c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 458c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid EigenSolver<MatrixType>::doComputeEigenvectors() 459c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 461c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index size = m_eivec.cols(); 462c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Scalar eps = NumTraits<Scalar>::epsilon(); 463c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 464c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // inefficient! this is already computed in RealSchur 465c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar norm(0); 466c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = 0; j < size; ++j) 467c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 468c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum(); 469c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 470c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 471c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Backsubstitute to find vectors of upper triangular form 4722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (norm == Scalar(0)) 473c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 474c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 475c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 476c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 477c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index n = size-1; n >= 0; n--) 478c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 479c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar p = m_eivalues.coeff(n).real(); 480c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar q = m_eivalues.coeff(n).imag(); 481c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 482c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Scalar vector 483c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (q == Scalar(0)) 484c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 485c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastr(0), lastw(0); 486c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n; 487c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matT.coeffRef(n,n) = Scalar(1); 489c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-1; i >= 0; i--) 490c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 491c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 492c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar r = m_matT.row(i).segment(l,n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 493c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 4942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_eivalues.coeff(i).imag() < Scalar(0)) 495c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 496c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 497c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastr = r; 498c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 499c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 500c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 501c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 5022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_eivalues.coeff(i).imag() == Scalar(0)) 503c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (w != Scalar(0)) 505c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / w; 506c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 507c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = -r / (eps * norm); 508c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 509c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else // Solve real equations 510c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 511c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 512c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 513c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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(); 514c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar t = (x * lastr - lastw * r) / denom; 515c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i,n) = t; 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(x) > abs(lastw)) 517c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-r - w * t) / x; 518c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 519c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-lastr - y * t) / lastw; 520c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 521c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 522c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar t = abs(m_matT.coeff(i,n)); 524c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 525c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.col(n).tail(size-i) /= t; 526c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 527c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 528c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 529c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else if (q < Scalar(0) && n > 0) // Complex vector 530c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 531c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar lastra(0), lastsa(0), lastw(0); 532c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index l = n-1; 533c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 534c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Last vector component imaginary so matrix is triangular 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(m_matT.coeff(n,n-1)) > abs(m_matT.coeff(n-1,n))) 536c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 537c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n-1) = q / m_matT.coeff(n,n-1); 538c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(n-1,n) = -(m_matT.coeff(n,n) - p) / m_matT.coeff(n,n-1); 539c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 540c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 541c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexScalar cc = ComplexScalar(Scalar(0),-m_matT.coeff(n-1,n)) / ComplexScalar(m_matT.coeff(n-1,n-1)-p,q); 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(n-1,n-1) = numext::real(cc); 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(n-1,n) = numext::imag(cc); 545c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 5462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matT.coeffRef(n,n-1) = Scalar(0); 5472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_matT.coeffRef(n,n) = Scalar(1); 548c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = n-2; i >= 0; i--) 549c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 550c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar ra = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n-1).segment(l, n-l+1)); 551c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar sa = m_matT.row(i).segment(l, n-l+1).dot(m_matT.col(n).segment(l, n-l+1)); 552c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar w = m_matT.coeff(i,i) - p; 553c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (m_eivalues.coeff(i).imag() < Scalar(0)) 555c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 556c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastw = w; 557c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastra = ra; 558c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath lastsa = sa; 559c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 560c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 561c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 562c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath l = i; 563c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (m_eivalues.coeff(i).imag() == RealScalar(0)) 564c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexScalar cc = ComplexScalar(-ra,-sa) / ComplexScalar(w,q); 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n-1) = numext::real(cc); 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n) = numext::imag(cc); 568c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 569c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 570c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 571c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Solve complex equations 572c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar x = m_matT.coeff(i,i+1); 573c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar y = m_matT.coeff(i+1,i); 574c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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; 575c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar vi = (m_eivalues.coeff(i).real() - p) * Scalar(2) * q; 5762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if ((vr == Scalar(0)) && (vi == Scalar(0))) 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez vr = eps * norm * (abs(w) + abs(q) + abs(x) + abs(y) + abs(lastw)); 578c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 5792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexScalar cc = ComplexScalar(x*lastra-lastw*ra+q*sa,x*lastsa-lastw*sa-q*ra) / ComplexScalar(vr,vi); 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n-1) = numext::real(cc); 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i,n) = numext::imag(cc); 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(x) > (abs(lastw) + abs(q))) 583c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 584c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n-1) = (-ra - w * m_matT.coeff(i,n-1) + q * m_matT.coeff(i,n)) / x; 585c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.coeffRef(i+1,n) = (-sa - w * m_matT.coeff(i,n) - q * m_matT.coeff(i,n-1)) / x; 586c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 587c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 588c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 5892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cc = ComplexScalar(-lastra-y*m_matT.coeff(i,n-1),-lastsa-y*m_matT.coeff(i,n)) / ComplexScalar(lastw,q); 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i+1,n-1) = numext::real(cc); 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_matT.coeffRef(i+1,n) = numext::imag(cc); 592c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 593c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 594c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 595c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Overflow control 5962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar t = numext::maxi<Scalar>(abs(m_matT.coeff(i,n-1)),abs(m_matT.coeff(i,n))); 597c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if ((eps * t) * t > Scalar(1)) 598c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matT.block(i, n-1, size-i, 2) /= t; 599c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 600c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 601c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 602c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 603c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // We handled a pair of complex conjugate eigenvalues, so need to skip them both 604c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath n--; 605c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 606c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath else 607c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 6082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(0 && "Internal bug in EigenSolver (INF or NaN has not been detected)"); // this should not happen 609c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 610c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 611c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 612c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Back transformation to get eigenvectors of original matrix 613c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index j = size-1; j >= 0; j--) 614c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 615c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_tmp.noalias() = m_eivec.leftCols(j+1) * m_matT.col(j).segment(0, j+1); 616c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(j) = m_tmp; 617c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 618c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 619c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 620c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 621c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 622c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_EIGENSOLVER_H 623