17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 42b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2012-2016 Gael Guennebaud <gael.guennebaud@inria.fr> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk> 62b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang// Copyright (C) 2016 Tobias Wood <tobias@spinicist.org.uk> 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_GENERALIZEDEIGENSOLVER_H 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_GENERALIZEDEIGENSOLVER_H 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#include "./RealQZ.h" 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez/** \eigenvalues_module \ingroup Eigenvalues_Module 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class GeneralizedEigenSolver 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Computes the generalized eigenvalues and eigenvectors of a pair of general matrices 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _MatrixType the type of the matrices of which we are computing the 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigen-decomposition; this is expected to be an instantiation of the Matrix 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * class template. Currently, only real matrices are supported. 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The generalized eigenvalues and eigenvectors of a matrix pair \f$ A \f$ and \f$ B \f$ are scalars 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f$ \lambda \f$ and vectors \f$ v \f$ such that \f$ Av = \lambda Bv \f$. If 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f$ D \f$ is a diagonal matrix with the eigenvalues on the diagonal, and 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f$ V \f$ is a matrix with the eigenvectors as its columns, then \f$ A V = 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * B V D \f$. The matrix \f$ V \f$ is almost always invertible, in which case we 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * have \f$ A = B V D V^{-1} \f$. This is called the generalized eigen-decomposition. 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The generalized eigenvalues and eigenvectors of a matrix pair may be complex, even when the 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrices are real. Moreover, the generalized eigenvalue might be infinite if the matrix B is 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * singular. To workaround this difficulty, the eigenvalues are provided as a pair of complex \f$ \alpha \f$ 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and real \f$ \beta \f$ such that: \f$ \lambda_i = \alpha_i / \beta_i \f$. If \f$ \beta_i \f$ is (nearly) zero, 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * then one can consider the well defined left eigenvalue \f$ \mu = \beta_i / \alpha_i\f$ such that: 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \f$ \mu_i A v_i = B v_i \f$, or even \f$ \mu_i u_i^T A = u_i^T B \f$ where \f$ u_i \f$ is 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * called the left eigenvector. 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Call the function compute() to compute the generalized eigenvalues and eigenvectors of 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * a given matrix pair. Alternatively, you can use the 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * GeneralizedEigenSolver(const MatrixType&, const MatrixType&, bool) constructor which computes the 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues and eigenvectors at construction time. Once the eigenvalue and 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvectors are computed, they can be retrieved with the eigenvalues() and 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvectors() functions. 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Here is an usage example of this class: 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include GeneralizedEigenSolver.cpp 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \verbinclude GeneralizedEigenSolver.out 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa MatrixBase::eigenvalues(), class ComplexEigenSolver, class SelfAdjointEigenSolver 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename _MatrixType> class GeneralizedEigenSolver 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Synonym for the template parameter \p _MatrixType. */ 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _MatrixType MatrixType; 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Options = MatrixType::Options, 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Scalar type for matrices of type #MatrixType. */ 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename NumTraits<Scalar>::Real RealScalar; 762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Complex scalar type for #MatrixType. 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is \c std::complex<Scalar> if #Scalar is real (e.g., 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \c float or \c double) and just \c Scalar if #Scalar is 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * complex. 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef std::complex<RealScalar> ComplexScalar; 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Type for vector of real scalar values eigenvalues as returned by betas(). 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is a column vector with entries of type #Scalar. 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The length of the vector is the size of #MatrixType. 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> VectorType; 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang /** \brief Type for vector of complex scalar values eigenvalues as returned by alphas(). 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is a column vector with entries of type #ComplexScalar. 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The length of the vector is the size of #MatrixType. 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ComplexVectorType; 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Expression type for the eigenvalues as returned by eigenvalues(). 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef CwiseBinaryOp<internal::scalar_quotient_op<ComplexScalar,Scalar>,ComplexVectorType,VectorType> EigenvalueType; 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Type for matrix of eigenvectors as returned by eigenvectors(). 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This is a square matrix with entries of type #ComplexScalar. 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The size is the same as the size of #MatrixType. 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime> EigenvectorsType; 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default constructor. 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The default constructor is useful in cases in which the user intends to 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * perform decompositions via EigenSolver::compute(const MatrixType&, bool). 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() for an example. 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1182b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang GeneralizedEigenSolver() 1192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang : m_eivec(), 1202b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_alphas(), 1212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_betas(), 1222b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_valuesOkay(false), 1232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_vectorsOkay(false), 1242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_realQZ() 1252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang {} 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default constructor with memory preallocation 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Like the default constructor but with preallocation of the internal data 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * according to the specified problem \a size. 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa GeneralizedEigenSolver() 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang explicit GeneralizedEigenSolver(Index size) 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_eivec(size, size), 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_alphas(size), 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_betas(size), 1372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_valuesOkay(false), 1382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_vectorsOkay(false), 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realQZ(size), 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_tmp(size) 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez {} 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Constructor; computes the generalized eigendecomposition of given matrix pair. 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Square matrix whose eigendecomposition is to be computed. 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Square matrix whose eigendecomposition is to be computed. 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] computeEigenvectors If true, both the eigenvectors and the 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues are computed; if false, only the eigenvalues are computed. 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This constructor calls compute() to compute the generalized eigenvalues 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and eigenvectors. 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez GeneralizedEigenSolver(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true) 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez : m_eivec(A.rows(), A.cols()), 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_alphas(A.cols()), 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_betas(A.cols()), 1592b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_valuesOkay(false), 1602b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_vectorsOkay(false), 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realQZ(A.cols()), 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_tmp(A.cols()) 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(A, B, computeEigenvectors); 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /* \brief Returns the computed generalized eigenvectors. 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * \returns %Matrix whose columns are the (possibly complex) right eigenvectors. 1702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang * i.e. the eigenvectors that solve (A - l*B)x = 0. The ordering matches the eigenvalues. 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre Either the constructor 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * GeneralizedEigenSolver(const MatrixType&,const MatrixType&, bool) or the member function 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * compute(const MatrixType&, const MatrixType& bool) has been called before, and 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \p computeEigenvectors was set to true (the default). 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa eigenvalues() 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EigenvectorsType eigenvectors() const { 1802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_vectorsOkay && "Eigenvectors for GeneralizedEigenSolver were not calculated."); 1812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang return m_eivec; 1822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns an expression of the computed generalized eigenvalues. 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns An expression of the column vector containing the eigenvalues. 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * It is a shortcut for \code this->alphas().cwiseQuotient(this->betas()); \endcode 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Not that betas might contain zeros. It is therefore not recommended to use this function, 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * but rather directly deal with the alphas and betas vectors. 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \pre Either the constructor 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * GeneralizedEigenSolver(const MatrixType&,const MatrixType&,bool) or the member function 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * compute(const MatrixType&,const MatrixType&,bool) has been called before. 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The eigenvalues are repeated according to their algebraic multiplicity, 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * so there are as many eigenvalues as rows in the matrix. The eigenvalues 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * are not sorted in any particular order. 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa alphas(), betas(), eigenvectors() 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez EigenvalueType eigenvalues() const 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return EigenvalueType(m_alphas,m_betas); 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns A const reference to the vectors containing the alpha values 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa betas(), eigenvalues() */ 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexVectorType alphas() const 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_alphas; 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \returns A const reference to the vectors containing the beta values 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This vector permits to reconstruct the j-th eigenvalues as alphas(i)/betas(j). 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa alphas(), eigenvalues() */ 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorType betas() const 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_valuesOkay && "GeneralizedEigenSolver is not initialized."); 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_betas; 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes generalized eigendecomposition of given matrix. 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Square matrix whose eigendecomposition is to be computed. 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Square matrix whose eigendecomposition is to be computed. 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] computeEigenvectors If true, both the eigenvectors and the 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * eigenvalues are computed; if false, only the eigenvalues are 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * computed. 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns Reference to \c *this 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This function computes the eigenvalues of the real matrix \p matrix. 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The eigenvalues() function can be used to retrieve them. If 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \p computeEigenvectors is true, then the eigenvectors are also computed 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and can be retrieved by calling eigenvectors(). 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The matrix is first reduced to real generalized Schur form using the RealQZ 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * class. The generalized Schur decomposition is then used to compute the eigenvalues 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * and eigenvectors. 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The cost of the computation is dominated by the cost of the 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * generalized Schur decomposition. 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This method reuses of the allocated data in the GeneralizedEigenSolver object. 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez GeneralizedEigenSolver& compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors = true); 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo info() const 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang eigen_assert(m_valuesOkay && "EigenSolver is not initialized."); 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_realQZ.info(); 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Sets the maximal number of iterations allowed. 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez GeneralizedEigenSolver& setMaxIterations(Index maxIters) 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realQZ.setMaxIterations(maxIters); 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez protected: 270a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 271a829215e078ace896f52702caa0c27608f40e3b0Miao Wang static void check_template_parameters() 272a829215e078ace896f52702caa0c27608f40e3b0Miao Wang { 273a829215e078ace896f52702caa0c27608f40e3b0Miao Wang EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 274a829215e078ace896f52702caa0c27608f40e3b0Miao Wang EIGEN_STATIC_ASSERT(!NumTraits<Scalar>::IsComplex, NUMERIC_TYPE_MUST_BE_REAL); 275a829215e078ace896f52702caa0c27608f40e3b0Miao Wang } 276a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 2772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang EigenvectorsType m_eivec; 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComplexVectorType m_alphas; 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez VectorType m_betas; 2802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang bool m_valuesOkay, m_vectorsOkay; 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ<MatrixType> m_realQZ; 2822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexVectorType m_tmp; 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez}; 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeztemplate<typename MatrixType> 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezGeneralizedEigenSolver<MatrixType>& 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos HernandezGeneralizedEigenSolver<MatrixType>::compute(const MatrixType& A, const MatrixType& B, bool computeEigenvectors) 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez{ 289a829215e078ace896f52702caa0c27608f40e3b0Miao Wang check_template_parameters(); 290a829215e078ace896f52702caa0c27608f40e3b0Miao Wang 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(A.cols() == A.rows() && B.cols() == A.rows() && B.cols() == B.rows()); 2942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Index size = A.cols(); 2952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_valuesOkay = false; 2962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_vectorsOkay = false; 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Reduce to generalized real Schur form: 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // A = Q S Z and B = Q T Z 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_realQZ.compute(A, B, computeEigenvectors); 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_realQZ.info() == Success) 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Resize storage 3032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_alphas.resize(size); 3042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_betas.resize(size); 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (computeEigenvectors) 3062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3072b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.resize(size,size); 3082b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_tmp.resize(size); 3092b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3102b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Aliases: 3122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Map<VectorType> v(reinterpret_cast<Scalar*>(m_tmp.data()), size); 3132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang ComplexVectorType &cv = m_tmp; 3142b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const MatrixType &mZ = m_realQZ.matrixZ(); 3152b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const MatrixType &mS = m_realQZ.matrixS(); 3162b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const MatrixType &mT = m_realQZ.matrixT(); 3172b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index i = 0; 3192b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang while (i < size) 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3212b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (i == size - 1 || mS.coeff(i+1, i) == Scalar(0)) 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3232b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Real eigenvalue 3242b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_alphas.coeffRef(i) = mS.diagonal().coeff(i); 3252b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_betas.coeffRef(i) = mT.diagonal().coeff(i); 3262b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (computeEigenvectors) 3272b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3282b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v.setConstant(Scalar(0.0)); 3292b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v.coeffRef(i) = Scalar(1.0); 3302b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // For singular eigenvalues do nothing more 3312b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if(abs(m_betas.coeffRef(i)) >= (std::numeric_limits<RealScalar>::min)()) 3322b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3332b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Non-singular eigenvalue 3342b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Scalar alpha = real(m_alphas.coeffRef(i)); 3352b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Scalar beta = m_betas.coeffRef(i); 3362b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (Index j = i-1; j >= 0; j--) 3372b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3382b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index st = j+1; 3392b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index sz = i-j; 3402b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) 3412b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3422b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // 2x2 block 3432b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<Scalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( v.segment(st,sz) ); 3442b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<Scalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1); 3452b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); 3462b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang j--; 3472b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3482b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang else 3492b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3502b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang v.coeffRef(j) = -v.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() / (beta*mS.coeffRef(j,j) - alpha*mT.coeffRef(j,j)); 3512b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3522b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3532b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3542b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i).real().noalias() = mZ.transpose() * v; 3552b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i).real().normalize(); 3562b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i).imag().setConstant(0); 3572b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ++i; 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3622b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // We need to extract the generalized eigenvalues of the pair of a general 2x2 block S and a positive diagonal 2x2 block T 3632b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Then taking beta=T_00*T_11, we can avoid any division, and alpha is the eigenvalues of A = (U^-1 * S * U) * diag(T_11,T_00): 3642b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3652b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // T = [a 0] 3662b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // [0 b] 3672b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang RealScalar a = mT.diagonal().coeff(i), 3682b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang b = mT.diagonal().coeff(i+1); 3692b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const RealScalar beta = m_betas.coeffRef(i) = m_betas.coeffRef(i+1) = a*b; 3702b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3712b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // ^^ NOTE: using diagonal()(i) instead of coeff(i,i) workarounds a MSVC bug. 3722b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<RealScalar,2,2> S2 = mS.template block<2,2>(i,i) * Matrix<Scalar,2,1>(b,a).asDiagonal(); 3732b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3742b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar p = Scalar(0.5) * (S2.coeff(0,0) - S2.coeff(1,1)); 3752b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Scalar z = sqrt(abs(p * p + S2.coeff(1,0) * S2.coeff(0,1))); 3762b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const ComplexScalar alpha = ComplexScalar(S2.coeff(1,1) + p, (beta > 0) ? z : -z); 3772b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_alphas.coeffRef(i) = conj(alpha); 3782b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_alphas.coeffRef(i+1) = alpha; 3792b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang 3802b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (computeEigenvectors) { 3812b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // Compute eigenvector in position (i+1) and then position (i) is just the conjugate 3822b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cv.setZero(); 3832b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cv.coeffRef(i+1) = Scalar(1.0); 3842b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // here, the "static_cast" workaound expression template issues. 3852b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cv.coeffRef(i) = -(static_cast<Scalar>(beta*mS.coeffRef(i,i+1)) - alpha*mT.coeffRef(i,i+1)) 3862b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang / (static_cast<Scalar>(beta*mS.coeffRef(i,i)) - alpha*mT.coeffRef(i,i)); 3872b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang for (Index j = i-1; j >= 0; j--) 3882b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3892b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index st = j+1; 3902b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang const Index sz = i+1-j; 3912b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang if (j > 0 && mS.coeff(j, j-1) != Scalar(0)) 3922b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang { 3932b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang // 2x2 block 3942b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<ComplexScalar, 2, 1> rhs = (alpha*mT.template block<2,Dynamic>(j-1,st,2,sz) - beta*mS.template block<2,Dynamic>(j-1,st,2,sz)) .lazyProduct( cv.segment(st,sz) ); 3952b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang Matrix<ComplexScalar, 2, 2> lhs = beta * mS.template block<2,2>(j-1,j-1) - alpha * mT.template block<2,2>(j-1,j-1); 3962b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cv.template segment<2>(j-1) = lhs.partialPivLu().solve(rhs); 3972b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang j--; 3982b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } else { 3992b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang cv.coeffRef(j) = cv.segment(st,sz).transpose().cwiseProduct(beta*mS.block(j,st,1,sz) - alpha*mT.block(j,st,1,sz)).sum() 4002b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang / (alpha*mT.coeffRef(j,j) - static_cast<Scalar>(beta*mS.coeffRef(j,j))); 4012b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4022b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4032b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i+1).noalias() = (mZ.transpose() * cv); 4042b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i+1).normalize(); 4052b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_eivec.col(i) = m_eivec.col(i+1).conjugate(); 4062b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez i += 2; 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4112b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_valuesOkay = true; 4122b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang m_vectorsOkay = computeEigenvectors; 4132b8756b6f1de65d3f8bffab45be6c44ceb7411fcMiao Wang } 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif // EIGEN_GENERALIZEDEIGENSOLVER_H 420