ComplexEigenSolver.h revision 7faaa9f3f0df9d23790277834d426c3d992ac3ba
165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// This file is part of Eigen, a lightweight C++ template library 265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// for linear algebra. 365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// 465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// Copyright (C) 2009 Claire Maurice 565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr> 665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// 865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// This Source Code Form is subject to the terms of the Mozilla 965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// Public License v. 2.0. If a copy of the MPL was not distributed 1065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 1165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 1265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich#ifndef EIGEN_COMPLEX_EIGEN_SOLVER_H 1365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich#define EIGEN_COMPLEX_EIGEN_SOLVER_H 1465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 1565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich#include "./ComplexSchur.h" 1665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 1765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichnamespace Eigen { 1865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 1965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich/** \eigenvalues_module \ingroup Eigenvalues_Module 2065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 2165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 2265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \class ComplexEigenSolver 2365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 2465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \brief Computes eigenvalues and eigenvectors of general complex matrices 2565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 2665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \tparam _MatrixType the type of the matrix of which we are 2765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * computing the eigendecomposition; this is expected to be an 2865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * instantiation of the Matrix class template. 2965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 3065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The eigenvalues and eigenvectors of a matrix \f$ A \f$ are scalars 3165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda v 3265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \f$. If \f$ D \f$ is a diagonal matrix with the eigenvalues on 3365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * the diagonal, and \f$ V \f$ is a matrix with the eigenvectors as 3465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * its columns, then \f$ A V = V D \f$. The matrix \f$ V \f$ is 3565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * almost always invertible, in which case we have \f$ A = V D V^{-1} 3665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \f$. This is called the eigendecomposition. 3765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 3865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The main function in this class is compute(), which computes the 3965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * eigenvalues and eigenvectors of a given function. The 4065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * documentation for that function contains an example showing the 4165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * main features of the class. 4265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 4365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \sa class EigenSolver, class SelfAdjointEigenSolver 4465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 4565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichtemplate<typename _MatrixType> class ComplexEigenSolver 4665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich{ 4765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich public: 4865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 4965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Synonym for the template parameter \p _MatrixType. */ 5065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef _MatrixType MatrixType; 5165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 5265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich enum { 5365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich RowsAtCompileTime = MatrixType::RowsAtCompileTime, 5465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ColsAtCompileTime = MatrixType::ColsAtCompileTime, 5565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich Options = MatrixType::Options, 5665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 5765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 5865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich }; 5965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 6065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Scalar type for matrices of type #MatrixType. */ 6165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef typename MatrixType::Scalar Scalar; 6265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef typename NumTraits<Scalar>::Real RealScalar; 6365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef typename MatrixType::Index Index; 6465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 6565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Complex scalar type for #MatrixType. 6665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 6765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This is \c std::complex<Scalar> if #Scalar is real (e.g., 6865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \c float or \c double) and just \c Scalar if #Scalar is 6965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * complex. 7065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 7165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef std::complex<RealScalar> ComplexScalar; 7265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 7365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Type for vector of eigenvalues as returned by eigenvalues(). 7465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 7565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This is a column vector with entries of type #ComplexScalar. 7665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The length of the vector is the size of #MatrixType. 7765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 7865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options&(~RowMajor), MaxColsAtCompileTime, 1> EigenvalueType; 7965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 8065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 8165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 8265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This is a square matrix with entries of type #ComplexScalar. 8365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The size is the same as the size of #MatrixType. 8465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 8565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorType; 8665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 8765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Default constructor. 8865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 8965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The default constructor is useful in cases in which the user intends to 9065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * perform decompositions via compute(). 9165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 9265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexEigenSolver() 9365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich : m_eivec(), 9465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivalues(), 9565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_schur(), 9665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_isInitialized(false), 9765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eigenvectorsOk(false), 9865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX() 9965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich {} 10065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 10165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Default Constructor with memory preallocation 10265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 10365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Like the default constructor but with preallocation of the internal data 10465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * according to the specified problem \a size. 10565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \sa ComplexEigenSolver() 10665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 10765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexEigenSolver(Index size) 10865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich : m_eivec(size, size), 10965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivalues(size), 11065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_schur(size), 11165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_isInitialized(false), 11265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eigenvectorsOk(false), 11365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX(size, size) 11465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich {} 11565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 11665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Constructor; computes eigendecomposition of given matrix. 11765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 11865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 11965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \param[in] computeEigenvectors If true, both the eigenvectors and the 12065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * eigenvalues are computed; if false, only the eigenvalues are 12165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * computed. 12265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 12365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This constructor calls compute() to compute the eigendecomposition. 12465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 12565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexEigenSolver(const MatrixType& matrix, bool computeEigenvectors = true) 12665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich : m_eivec(matrix.rows(),matrix.cols()), 12765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivalues(matrix.cols()), 12865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_schur(matrix.rows()), 12965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_isInitialized(false), 13065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eigenvectorsOk(false), 13165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX(matrix.rows(),matrix.cols()) 13265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 13365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich compute(matrix, computeEigenvectors); 13465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 13565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 13665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Returns the eigenvectors of given matrix. 13765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 13865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \returns A const reference to the matrix whose columns are the eigenvectors. 13965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 14065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \pre Either the constructor 14165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 14265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * function compute(const MatrixType& matrix, bool) has been called before 14365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * to compute the eigendecomposition of a matrix, and 14465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \p computeEigenvectors was set to true (the default). 14565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 14665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This function returns a matrix whose columns are the eigenvectors. Column 14765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \f$ k \f$ is an eigenvector corresponding to eigenvalue number \f$ k 14865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \f$ as returned by eigenvalues(). The eigenvectors are normalized to 14965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * have (Euclidean) norm equal to one. The matrix returned by this 15065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * function is the matrix \f$ V \f$ in the eigendecomposition \f$ A = V D 15165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * V^{-1} \f$, if it exists. 15265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 15365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Example: \include ComplexEigenSolver_eigenvectors.cpp 15465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Output: \verbinclude ComplexEigenSolver_eigenvectors.out 15565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 15665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich const EigenvectorType& eigenvectors() const 15765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 15865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 15965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich eigen_assert(m_eigenvectorsOk && "The eigenvectors have not been computed together with the eigenvalues."); 16065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return m_eivec; 16165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 16265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 16365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Returns the eigenvalues of given matrix. 16465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 16565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \returns A const reference to the column vector containing the eigenvalues. 16665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 16765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \pre Either the constructor 16865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * ComplexEigenSolver(const MatrixType& matrix, bool) or the member 16965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * function compute(const MatrixType& matrix, bool) has been called before 17065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * to compute the eigendecomposition of a matrix. 17165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 17265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This function returns a column vector containing the 17365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * eigenvalues. Eigenvalues are repeated according to their 17465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * algebraic multiplicity, so there are as many eigenvalues as 17565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * rows in the matrix. The eigenvalues are not sorted in any particular 17665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * order. 17765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 17865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Example: \include ComplexEigenSolver_eigenvalues.cpp 17965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Output: \verbinclude ComplexEigenSolver_eigenvalues.out 18065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 18165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich const EigenvalueType& eigenvalues() const 18265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 18365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 18465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return m_eivalues; 18565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 18665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 18765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Computes eigendecomposition of given matrix. 18865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 18965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \param[in] matrix Square matrix whose eigendecomposition is to be computed. 19065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \param[in] computeEigenvectors If true, both the eigenvectors and the 19165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * eigenvalues are computed; if false, only the eigenvalues are 19265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * computed. 19365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \returns Reference to \c *this 19465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 19565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * This function computes the eigenvalues of the complex matrix \p matrix. 19665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The eigenvalues() function can be used to retrieve them. If 19765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \p computeEigenvectors is true, then the eigenvectors are also computed 19865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * and can be retrieved by calling eigenvectors(). 19965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 20065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The matrix is first reduced to Schur form using the 20165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * ComplexSchur class. The Schur decomposition is then used to 20265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * compute the eigenvalues and eigenvectors. 20365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 20465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * The cost of the computation is dominated by the cost of the 20565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Schur decomposition, which is \f$ O(n^3) \f$ where \f$ n \f$ 20665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * is the size of the matrix. 20765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 20865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Example: \include ComplexEigenSolver_compute.cpp 20965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * Output: \verbinclude ComplexEigenSolver_compute.out 21065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 21165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexEigenSolver& compute(const MatrixType& matrix, bool computeEigenvectors = true); 21265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 21365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Reports whether previous computation was successful. 21465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * 21565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 21665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich */ 21765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComputationInfo info() const 21865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 21965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich eigen_assert(m_isInitialized && "ComplexEigenSolver is not initialized."); 22065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return m_schur.info(); 22165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 22265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 22365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Sets the maximum number of iterations allowed. */ 22465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexEigenSolver& setMaxIterations(Index maxIters) 22565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 22665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_schur.setMaxIterations(maxIters); 22765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return *this; 22865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 22965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 23065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich /** \brief Returns the maximum number of iterations. */ 23165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich Index getMaxIterations() 23265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 23365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return m_schur.getMaxIterations(); 23465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 23565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 23665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich protected: 23765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich EigenvectorType m_eivec; 23865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich EigenvalueType m_eivalues; 23965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexSchur<MatrixType> m_schur; 24065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich bool m_isInitialized; 24165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich bool m_eigenvectorsOk; 24265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich EigenvectorType m_matX; 24365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 24465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich private: 24565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich void doComputeEigenvectors(const RealScalar& matrixnorm); 24665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich void sortEigenvalues(bool computeEigenvectors); 24765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich}; 24865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 24965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 25065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichtemplate<typename MatrixType> 25165de34233da93a3d65c00b8aad3ff9aad44c57deNick KralevichComplexEigenSolver<MatrixType>& 25265de34233da93a3d65c00b8aad3ff9aad44c57deNick KralevichComplexEigenSolver<MatrixType>::compute(const MatrixType& matrix, bool computeEigenvectors) 25365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich{ 25465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // this code is inspired from Jampack 25565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich eigen_assert(matrix.cols() == matrix.rows()); 25665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 25765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // Do a complex Schur decomposition, A = U T U^* 25865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // The eigenvalues are on the diagonal of T. 25965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_schur.compute(matrix, computeEigenvectors); 26065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 26165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if(m_schur.info() == Success) 26265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 26365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivalues = m_schur.matrixT().diagonal(); 26465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if(computeEigenvectors) 26565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich doComputeEigenvectors(matrix.norm()); 26665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich sortEigenvalues(computeEigenvectors); 26765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 26865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 26965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_isInitialized = true; 27065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eigenvectorsOk = computeEigenvectors; 27165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich return *this; 27265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich} 27365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 27465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 27565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichtemplate<typename MatrixType> 27665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichvoid ComplexEigenSolver<MatrixType>::doComputeEigenvectors(const RealScalar& matrixnorm) 27765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich{ 27865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich const Index n = m_eivalues.size(); 27965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 28065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // Compute X such that T = X D X^(-1), where D is the diagonal of T. 28165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // The matrix X is unit triangular. 28265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX = EigenvectorType::Zero(n, n); 28365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich for(Index k=n-1 ; k>=0 ; k--) 28465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 28565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX.coeffRef(k,k) = ComplexScalar(1.0,0.0); 28665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // Compute X(i,k) using the (i,k) entry of the equation X T = D X 28765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich for(Index i=k-1 ; i>=0 ; i--) 28865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 28965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX.coeffRef(i,k) = -m_schur.matrixT().coeff(i,k); 29065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if(k-i-1>0) 29165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 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(); 29265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich ComplexScalar z = m_schur.matrixT().coeff(i,i) - m_schur.matrixT().coeff(k,k); 29365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if(z==ComplexScalar(0)) 29465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 29565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // If the i-th and k-th eigenvalue are equal, then z equals 0. 29665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // Use a small value instead, to prevent division by zero. 29765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich numext::real_ref(z) = NumTraits<RealScalar>::epsilon() * matrixnorm; 29865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 29965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_matX.coeffRef(i,k) = m_matX.coeff(i,k) / z; 30065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 30165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 30265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 30365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // Compute V as V = U X; now A = U T U^* = U X D X^(-1) U^* = V D V^(-1) 30465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivec.noalias() = m_schur.matrixU() * m_matX; 30565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich // .. and normalize the eigenvectors 30665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich for(Index k=0 ; k<n ; k++) 30765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 30865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivec.col(k).normalize(); 30965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 31065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich} 31165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 31265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 31365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichtemplate<typename MatrixType> 31465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevichvoid ComplexEigenSolver<MatrixType>::sortEigenvalues(bool computeEigenvectors) 31565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich{ 31665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich const Index n = m_eivalues.size(); 31765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich for (Index i=0; i<n; i++) 31865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 31965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich Index k; 32065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivalues.cwiseAbs().tail(n-i).minCoeff(&k); 32165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if (k != 0) 32265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich { 32365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich k += i; 32465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich std::swap(m_eivalues[k],m_eivalues[i]); 32565de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich if(computeEigenvectors) 32665de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich m_eivec.col(i).swap(m_eivec.col(k)); 32765de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 32865de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich } 32965de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich} 33065de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 33165de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich} // end namespace Eigen 33265de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich 33365de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich#endif // EIGEN_COMPLEX_EIGEN_SOLVER_H 33465de34233da93a3d65c00b8aad3ff9aad44c57deNick Kralevich