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> 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 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; 632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 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 */ 1072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit 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 */ 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 1262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit ComplexEigenSolver(const EigenBase<InputType>& matrix, bool computeEigenvectors = true) 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_eivec(matrix.rows(),matrix.cols()), 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues(matrix.cols()), 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_schur(matrix.rows()), 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false), 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk(false), 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX(matrix.rows(),matrix.cols()) 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 1342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang compute(matrix.derived(), computeEigenvectors); 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvectors of given matrix. 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the matrix whose columns are the eigenvectors. 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function compute(const MatrixType& matrix, bool) has been called before 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigendecomposition of a matrix, and 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors was set to true (the default). 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a matrix whose columns are the eigenvectors. Column 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ as returned by eigenvalues(). The eigenvectors are normalized to 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * have (Euclidean) norm equal to one. The matrix returned by this 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * V^{-1} \f$, if it exists. 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_eigenvectors.cpp 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_eigenvectors.out 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvectorType& eigenvectors() const 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivec; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the eigenvalues of given matrix. 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns A const reference to the column vector containing the eigenvalues. 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * function compute(const MatrixType& matrix, bool) has been called before 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * to compute the eigendecomposition of a matrix. 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a column vector containing the 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues. Eigenvalues are repeated according to their 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * algebraic multiplicity, so there are as many eigenvalues as 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * rows in the matrix. The eigenvalues are not sorted in any particular 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * order. 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_eigenvalues.cpp 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_eigenvalues.out 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const EigenvalueType& eigenvalues() const 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_eivalues; 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes eigendecomposition of given matrix. 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] computeEigenvectors If true, both the eigenvectors and the 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * eigenvalues are computed; if false, only the eigenvalues are 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed. 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function computes the eigenvalues of the complex matrix \p matrix. 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The eigenvalues() function can be used to retrieve them. If 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \p computeEigenvectors is true, then the eigenvectors are also computed 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * and can be retrieved by calling eigenvectors(). 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The matrix is first reduced to Schur form using the 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * ComplexSchur class. The Schur decomposition is then used to 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * compute the eigenvalues and eigenvectors. 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The cost of the computation is dominated by the cost of the 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * is the size of the matrix. 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include ComplexEigenSolver_compute.cpp 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude ComplexEigenSolver_compute.out 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 2122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang template<typename InputType> 2132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexEigenSolver& compute(const EigenBase<InputType>& matrix, bool computeEigenvectors = true); 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reports whether previous computation was successful. 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComputationInfo info() const 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_schur.info(); 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Sets the maximum number of iterations allowed. */ 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexEigenSolver& setMaxIterations(Index maxIters) 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_schur.setMaxIterations(maxIters); 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns the maximum number of iterations. */ 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index getMaxIterations() 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_schur.getMaxIterations(); 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 239a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 240a829215e078ace896f52702caa0c27608f40e3b0Miao Wang static void check_template_parameters() 241a829215e078ace896f52702caa0c27608f40e3b0Miao Wang { 242a829215e078ace896f52702caa0c27608f40e3b0Miao Wang EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 243a829215e078ace896f52702caa0c27608f40e3b0Miao Wang } 244a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorType m_eivec; 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvalueType m_eivalues; 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexSchur<MatrixType> m_schur; 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_eigenvectorsOk; 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath EigenvectorType m_matX; 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 2532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang void doComputeEigenvectors(RealScalar matrixnorm); 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath void sortEigenvalues(bool computeEigenvectors); 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 2592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangtemplate<typename InputType> 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezComplexEigenSolver<MatrixType>& 2612b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao WangComplexEigenSolver<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeEigenvectors) 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 263a829215e078ace896f52702caa0c27608f40e3b0Miao Wang check_template_parameters(); 264a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // this code is inspired from Jampack 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(matrix.cols() == matrix.rows()); 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Do a complex Schur decomposition, A = U T U^* 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The eigenvalues are on the diagonal of T. 2702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_schur.compute(matrix.derived(), computeEigenvectors); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(m_schur.info() == Success) 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues = m_schur.matrixT().diagonal(); 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigenvectors) 2762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang doComputeEigenvectors(m_schur.matrixT().norm()); 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath sortEigenvalues(computeEigenvectors); 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eigenvectorsOk = computeEigenvectors; 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 2872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wangvoid ComplexEigenSolver<MatrixType>::doComputeEigenvectors(RealScalar matrixnorm) 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index n = m_eivalues.size(); 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 2912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang matrixnorm = numext::maxi(matrixnorm,(std::numeric_limits<RealScalar>::min)()); 2922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute X such that T = X D X^(-1), where D is the diagonal of T. 294c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // The matrix X is unit triangular. 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX = EigenvectorType::Zero(n, n); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k=n-1 ; k>=0 ; k--) 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute X(i,k) using the (i,k) entry of the equation X T = D X 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index i=k-1 ; i>=0 ; i--) 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(k-i-1>0) 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan 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(); 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(z==ComplexScalar(0)) 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // If the i-th and k-th eigenvalue are equal, then z equals 0. 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Use a small value instead, to prevent division by zero. 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 316c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.noalias() = m_schur.matrixU() * m_matX; 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // .. and normalize the eigenvectors 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for(Index k=0 ; k<n ; k++) 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(k).normalize(); 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const Index n = m_eivalues.size(); 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i=0; i<n; i++) 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index k; 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (k != 0) 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath k += i; 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath std::swap(m_eivalues[k],m_eivalues[i]); 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(computeEigenvectors) 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_eivec.col(i).swap(m_eivec.col(k)); 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H 347