1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Claire Maurice 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_COMPLEX_EIGEN_SOLVER_H 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#include "./ComplexSchur.h" 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class ComplexEigenSolver 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Computes eigenvalues and eigenvectors of general complex matrices 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computing the eigendecomposition; this is expected to be an 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * instantiation of the Matrix class template. 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$. If \f$ D \f$ is a diagonal matrix with the eigenvalues on 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * almost always invertible, in which case we have \f$ A = V D V^{-1} 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$. This is called the eigendecomposition. 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The main function in this class is compute(), which computes the 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues and eigenvectors of a given function. The 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * documentation for that function contains an example showing the 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * main features of the class. 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class EigenSolver, class SelfAdjointEigenSolver 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class ComplexEigenSolver 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Synonym for the template parameter \p _MatrixType. */ 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RowsAtCompileTime = MatrixType::RowsAtCompileTime, 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ColsAtCompileTime = MatrixType::ColsAtCompileTime, 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Scalar type for matrices of type #MatrixType. */ 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Complex scalar type for #MatrixType. 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is \c std::complex<Scalar> if #Scalar is real (e.g., 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \c float or \c double) and just \c Scalar if #Scalar is 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * complex. 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef std::complex<RealScalar> ComplexScalar; 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a column vector with entries of type #ComplexScalar. 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The length of the vector is the size of #MatrixType. 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType; 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is a square matrix with entries of type #ComplexScalar. 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The size is the same as the size of #MatrixType. 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 85c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor. 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via compute(). 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexEigenSolver() 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(), 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(), 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_schur(), 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX() 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default Constructor with memory preallocation 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Like the default constructor but with preallocation of the internal data 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * according to the specified problem \a size. 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa ComplexEigenSolver() 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexEigenSolver(Index size) 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(size, size), 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(size), 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_schur(size), 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX(size, size) 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath {} 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes eigendecomposition of given matrix. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute() to compute the eigendecomposition. 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(matrix.rows(),matrix.cols()), 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(matrix.cols()), 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_schur(matrix.rows()), 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX(matrix.rows(),matrix.cols()) 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath compute(matrix, computeEigenvectors); 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvectors of given matrix. 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix whose columns are the eigenvectors. 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function compute(const MatrixType& matrix, bool) has been called before 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigendecomposition of a matrix, and 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a matrix whose columns are the eigenvectors. Column 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ as returned by eigenvalues(). The eigenvectors are normalized to 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * have (Euclidean) norm equal to one. The matrix returned by this 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * V^{-1} \f$, if it exists. 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_eigenvectors.cpp 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_eigenvectors.out 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvectorType& eigenvectors() const 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivec; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvalues of given matrix. 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the column vector containing the eigenvalues. 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function compute(const MatrixType& matrix, bool) has been called before 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigendecomposition of a matrix. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a column vector containing the 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues. Eigenvalues are repeated according to their 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * algebraic multiplicity, so there are as many eigenvalues as 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * rows in the matrix. The eigenvalues are not sorted in any particular 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * order. 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_eigenvalues.cpp 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_eigenvalues.out 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvalueType& eigenvalues() const 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivalues; 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes eigendecomposition of given matrix. 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function computes the eigenvalues of the complex matrix \p matrix. 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve them. If 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors is true, then the eigenvectors are also computed 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and can be retrieved by calling eigenvectors(). 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix is first reduced to Schur form using the 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur class. The Schur decomposition is then used to 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute the eigenvalues and eigenvectors. 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The cost of the computation is dominated by the cost of the 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is the size of the matrix. 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_compute.cpp 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_compute.out 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_schur.info(); 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorType m_eivec; 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvalueType m_eivalues; 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexSchur<MatrixType> m_schur; 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_eigenvectorsOk; 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorType m_matX; 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void doComputeEigenvectors(RealScalar matrixnorm); 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void sortEigenvalues(bool computeEigenvectors); 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan KamathComplexEigenSolver<MatrixType>& ComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // this code is inspired from Jampack 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath assert(matrix.cols() == matrix.rows()); 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Do a complex Schur decomposition, A = U T U^* 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The eigenvalues are on the diagonal of T. 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_schur.compute(matrix, computeEigenvectors); 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_schur.info() == Success) 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues = m_schur.matrixT().diagonal(); 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigenvectors) 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath doComputeEigenvectors(matrix.norm()); 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sortEigenvalues(computeEigenvectors); 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk = computeEigenvectors; 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index n = m_eivalues.size(); 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute X such that T = X D X^(-1), where D is the diagonal of T. 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The matrix X is unit triangular. 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX = EigenvectorType::Zero(n, n); 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k=n-1 ; k>=0 ; k--) 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute X(i,k) using the (i,k) entry of the equation X T = D X 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=k-1 ; i>=0 ; i--) 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k-i-1>0) 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(i,k) -= (m_schur.matrixT().row(i).segment(i+1,k-i-1) * m_matX.col(k).segment(i+1,k-i-1)).value(); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(z==ComplexScalar(0)) 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // If the i-th and k-th eigenvalue are equal, then z equals 0. 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Use a small value instead, to prevent division by zero. 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath internal::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.noalias() = m_schur.matrixU() * m_matX; 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // .. and normalize the eigenvectors 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k=0 ; k<n ; k++) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(k).normalize(); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index n = m_eivalues.size(); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; i++) 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != 0) 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k += i; 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(m_eivalues[k],m_eivalues[i]); 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigenvectors) 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(i).swap(m_eivec.col(k)); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H 320