17faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This file is part of Eigen, a lightweight C++ template library 27faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// for linear algebra. 37faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 47faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru> 57faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// 67faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// This Source Code Form is subject to the terms of the Mozilla 77faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// Public License v. 2.0. If a copy of the MPL was not distributed 87faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 97faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#ifndef EIGEN_REAL_QZ_H 117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#define EIGEN_REAL_QZ_H 127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandeznamespace Eigen { 147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \eigenvalues_module \ingroup Eigenvalues_Module 167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \class RealQZ 197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \brief Performs a real QZ decomposition of a pair of square matrices 217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \tparam _MatrixType the type of the matrix of which we are computing the 237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * real QZ decomposition; this is expected to be an instantiation of the 247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Matrix class template. 257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Given a real square matrices A and B, this class computes the real QZ 277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are 287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * real orthogonal matrixes, T is upper-triangular matrix, and S is upper 297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * quasi-triangular matrix. An orthogonal matrix is a matrix whose 307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular 317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrix is a block-triangular matrix whose diagonal consists of 1-by-1 327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * blocks and 2-by-2 blocks where further reduction is impossible due to 337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * complex eigenvalues. 347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from 367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1x1 and 2x2 blocks on the diagonals of S and T. 377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Call the function compute() to compute the real QZ decomposition of a 397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * given pair of matrices. Alternatively, you can use the 407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ) 417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * constructor which computes the real QZ decomposition at construction 427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * time. Once the decomposition is computed, you can use the matrixS(), 437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices 447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * S, T, Q and Z in the decomposition. If computeQZ==false, some time 457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * is saved by not computing matrices Q and Z. 467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Example: \include RealQZ_compute.cpp 487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * Output: \include RealQZ_compute.out 497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \note The implementation is based on the algorithm in "Matrix Computations" 517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for 527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart. 537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver 557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename _MatrixType> class RealQZ 587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez public: 607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef _MatrixType MatrixType; 617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez enum { 627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RowsAtCompileTime = MatrixType::RowsAtCompileTime, 637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ColsAtCompileTime = MatrixType::ColsAtCompileTime, 647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Options = MatrixType::Options, 657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; 687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Scalar Scalar; 697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef typename MatrixType::Index Index; 717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Default constructor. 767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param [in] size Positive integer, size of the matrix whose QZ decomposition will be computed. 787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * The default constructor is useful in cases in which the user intends to 807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * perform decompositions via compute(). The \p size parameter is only 817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * used as a hint. It is not an error to give a wrong \p size, but it may 827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * impair performance. 837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \sa compute() for an example. 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : 877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S(size, size), 887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T(size, size), 897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q(size, size), 907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z(size, size), 917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workspace(size*2), 927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters(400), 937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false) 947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { } 957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Constructor; computes real QZ decomposition of given matrices 977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Matrix A. 997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Matrix B. 1007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] computeQZ If false, A and Z are not computed. 1017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * This constructor calls compute() to compute the QZ decomposition. 1037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) : 1057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S(A.rows(),A.cols()), 1067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T(A.rows(),A.cols()), 1077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q(A.rows(),A.cols()), 1087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z(A.rows(),A.cols()), 1097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workspace(A.rows()*2), 1107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters(400), 1117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized(false) { 1127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez compute(A, B, computeQZ); 1137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns matrix Q in the QZ decomposition. 1167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the matrix Q. 1187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixType& matrixQ() const { 1207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); 1227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_Q; 1237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns matrix Z in the QZ decomposition. 1267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the matrix Z. 1287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixType& matrixZ() const { 1307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); 1327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_Z; 1337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns matrix S in the QZ decomposition. 1367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the matrix S. 1387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixType& matrixS() const { 1407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_S; 1427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns matrix S in the QZ decomposition. 1457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns A const reference to the matrix S. 1477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const MatrixType& matrixT() const { 1497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_T; 1517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Computes QZ decomposition of given matrix. 1547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] A Matrix A. 1567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] B Matrix B. 1577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \param[in] computeQZ If false, A and Z are not computed. 1587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns Reference to \c *this 1597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true); 1617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Reports whether previous computation was successful. 1637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * 1647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 1657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo info() const 1677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_info; 1707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \brief Returns number of performed QR-like iterations. 1737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index iterations() const 1757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(m_isInitialized && "RealQZ is not initialized."); 1777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return m_global_iter; 1787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** Sets the maximal number of iterations allowed to converge to one eigenvalue 1817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez * or decouple the problem. 1827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez */ 1837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ& setMaxIterations(Index maxIters) 1847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 1857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_maxIters = maxIters; 1867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 1877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 1887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez private: 1907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 1917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez MatrixType m_S, m_T, m_Q, m_Z; 1927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix<Scalar,Dynamic,1> m_workspace; 1937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ComputationInfo m_info; 1947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_maxIters; 1957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_isInitialized; 1967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez bool m_computeQZ; 1977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar m_normOfT, m_normOfS; 1987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index m_global_iter; 1997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,3,1> Vector3s; 2017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,2,1> Vector2s; 2027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef Matrix<Scalar,2,2> Matrix2s; 2037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef JacobiRotation<Scalar> JRs; 2047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void hessenbergTriangular(); 2067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void computeNorms(); 2077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index findSmallSubdiagEntry(Index iu); 2087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index findSmallDiagEntry(Index f, Index l); 2097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void splitOffTwoRows(Index i); 2107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void pushDownZero(Index z, Index f, Index l); 2117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void step(Index f, Index l, Index iter); 2127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez }; // RealQZ 2147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal Reduces S and T to upper Hessenberg - triangular form */ 2167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 2177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez void RealQZ<MatrixType>::hessenbergTriangular() 2187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index dim = m_S.cols(); 2217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // perform QR decomposition of T, overwrite T with R, save Q 2237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez HouseholderQR<MatrixType> qrT(m_T); 2247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T = qrT.matrixQR(); 2257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.template triangularView<StrictlyLower>().setZero(); 2267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q = qrT.householderQ(); 2277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // overwrite S with Q* S 2287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheLeft(m_Q.adjoint()); 2297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // init Z as Identity 2307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 2317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z = MatrixType::Identity(dim,dim); 2327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // reduce S to upper Hessenberg with Givens rotations 2337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j=0; j<=dim-3; j++) { 2347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index i=dim-1; i>=j+2; i--) { 2357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JRs G; 2367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // kill S(i,j) 2377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(m_S.coeff(i,j) != 0) 2387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j)); 2407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(i,j) = Scalar(0.0); 2417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint()); 2427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint()); 2437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Q 2457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 2467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q.applyOnTheRight(i-1,i,G); 2477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // kill T(i,i-1) 2487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if(m_T.coeff(i,i-1)!=Scalar(0)) 2497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i)); 2517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(i,i-1) = Scalar(0.0); 2527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheRight(i,i-1,G); 2537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.topRows(i).applyOnTheRight(i,i-1,G); 2547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Z 2567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 2577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(i,i-1,G.adjoint()); 2587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */ 2637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 2647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline void RealQZ<MatrixType>::computeNorms() 2657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index size = m_S.cols(); 2677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_normOfS = Scalar(0.0); 2687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_normOfT = Scalar(0.0); 2697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index j = 0; j < size; ++j) 2707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); 2727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum(); 2737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ 2787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 2797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) 2807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 2827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index res = iu; 2837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while (res > 0) 2847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res)); 2867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (s == Scalar(0.0)) 2877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez s = m_normOfS; 2887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s) 2897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 2907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res--; 2917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return res; 2937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 2957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */ 2967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 2977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) 2987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 2997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 3007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index res = l; 3017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while (res >= f) { 3027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT) 3037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez break; 3047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez res--; 3057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return res; 3077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */ 3107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 3117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) 3127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 3147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::sqrt; 3157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index dim=m_S.cols(); 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (abs(m_S.coeff(i+1,i)==Scalar(0))) 3177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return; 3187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index z = findSmallDiagEntry(i,i+1); 3197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (z==i-1) 3207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // block of (S T^{-1}) 3227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>(). 3237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template solve<OnTheRight>(m_S.template block<2,2>(i,i)); 3247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1)); 3257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar q = p*p + STi(1,0)*STi(0,1); 3267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (q>=0) { 3277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar z = sqrt(q); 3287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // one QR-like iteration for ABi - lambda I 3297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // is enough - when we know exact eigenvalue in advance, 3307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // convergence is immediate 3317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JRs G; 3327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (p>=0) 3337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(p + z, STi(1,0)); 3347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(p - z, STi(1,0)); 3367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); 3377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); 3387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Q 3397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 3407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q.applyOnTheRight(i,i+1,G); 3417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i)); 3437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.topRows(i+2).applyOnTheRight(i+1,i,G); 3447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.topRows(i+2).applyOnTheRight(i+1,i,G); 3457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Z 3467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 3477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(i+1,i,G.adjoint()); 3487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(i+1,i) = Scalar(0.0); 3507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(i+1,i) = Scalar(0.0); 3517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 3547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez pushDownZero(z,i,i+1); 3567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */ 3607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 3617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l) 3627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JRs G; 3647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index dim = m_S.cols(); 3657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index zz=z; zz<l; zz++) 3667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // push 0 down 3687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index firstColS = zz>f ? (zz-1) : zz; 3697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1)); 3707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint()); 3717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint()); 3727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(zz+1,zz+1) = Scalar(0.0); 3737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Q 3747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 3757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q.applyOnTheRight(zz,zz+1,G); 3767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // kill S(zz+1, zz-1) 3777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (zz>f) 3787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 3797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1)); 3807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G); 3817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G); 3827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(zz+1,zz-1) = Scalar(0.0); 3837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Z 3847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 3857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(zz,zz-1,G.adjoint()); 3867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // finally kill S(l,l-1) 3897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1)); 3907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheRight(l,l-1,G); 3917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.applyOnTheRight(l,l-1,G); 3927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(l,l-1)=Scalar(0.0); 3937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Z 3947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 3957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(l,l-1,G.adjoint()); 3967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 3977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 3987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez /** \internal QR-like iterative step for block f..l */ 3997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 4007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) 4017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez using std::abs; 4037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index dim = m_S.cols(); 4047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // x, y, z 4067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar x, y, z; 4077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (iter==10) 4087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Wilkinson ad hoc shift 4107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Scalar 4117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1), 4127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1), 4137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b12=m_T.coeff(f+0,f+1), 4147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b11i=Scalar(1.0)/m_T.coeff(f+0,f+0), 4157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b22i=Scalar(1.0)/m_T.coeff(f+1,f+1), 4167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a87=m_S.coeff(l-1,l-2), 4177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a98=m_S.coeff(l-0,l-1), 4187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b77i=Scalar(1.0)/m_T.coeff(l-2,l-2), 4197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b88i=Scalar(1.0)/m_T.coeff(l-1,l-1); 4207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar ss = abs(a87*b77i) + abs(a98*b88i), 4217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez lpl = Scalar(1.5)*ss, 4227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez ll = ss*ss; 4237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i 4247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez - a11*a21*b12*b11i*b11i*b22i; 4257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i 4267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez - a21*a21*b12*b11i*b11i*b22i; 4277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = a21*a32*b11i*b22i; 4287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (iter==16) 4307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // another exceptional shift 4327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = m_S.coeff(f,f)/m_T.coeff(f,f)-m_S.coeff(l,l)/m_T.coeff(l,l) + m_S.coeff(l,l-1)*m_T.coeff(l-1,l) / 4337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez (m_T.coeff(l-1,l-1)*m_T.coeff(l,l)); 4347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = m_S.coeff(f+1,f)/m_T.coeff(f,f); 4357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = 0; 4367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (iter>23 && !(iter%8)) 4387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // extremely exceptional shift 4407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = internal::random<Scalar>(-1.0,1.0); 4417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = internal::random<Scalar>(-1.0,1.0); 4427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = internal::random<Scalar>(-1.0,1.0); 4437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 4457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1 4477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where 4487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // U and V are 2x2 bottom right sub matrices of A and B. Thus: 4497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1) 4507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // = AB^-1AB^-1 + det(M) - tr(M)(AB^-1) 4517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Since we are only interested in having x, y, z with a correct ratio, we have: 4527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Scalar 4537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a11 = m_S.coeff(f,f), a12 = m_S.coeff(f,f+1), 4547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a21 = m_S.coeff(f+1,f), a22 = m_S.coeff(f+1,f+1), 4557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a32 = m_S.coeff(f+2,f+1), 4567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l), 4587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez a98 = m_S.coeff(l,l-1), a99 = m_S.coeff(l,l), 4597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b11 = m_T.coeff(f,f), b12 = m_T.coeff(f,f+1), 4617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b22 = m_T.coeff(f+1,f+1), 4627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l), 4647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez b99 = m_T.coeff(l,l); 4657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21) 4677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez + a12/b22 - (a11/b11)*(b12/b22); 4687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99); 4697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = a32/b22; 4707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 4717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez JRs G; 4737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez for (Index k=f; k<=l-2; k++) 4757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // variables for Householder reflections 4777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Vector2s essential2; 4787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Scalar tau, beta; 4797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Vector3s hr(x,y,z); 4817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1) 4837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hr.makeHouseholderInPlace(tau, beta); 4847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez essential2 = hr.template bottomRows<2>(); 4857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index fc=(std::max)(k-1,Index(0)); // first col to update 4867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); 4877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); 4887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 4897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data()); 4907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (k>f) 4917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0); 4927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 4937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k) 4947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1); 4957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hr.makeHouseholderInPlace(tau, beta); 4967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez essential2 = hr.template bottomRows<2>(); 4977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 4987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index lr = (std::min)(k+4,dim); // last row to update 4997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr); 5007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // S 5017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2; 5027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp += m_S.col(k+2).head(lr); 5037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.col(k+2).head(lr) -= tau*tmp; 5047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); 5057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // T 5067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2; 5077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp += m_T.col(k+2).head(lr); 5087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.col(k+2).head(lr) -= tau*tmp; 5097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); 5107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 5127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Z 5147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim); 5157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k)); 5167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez tmp += m_Z.row(k+2); 5177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.row(k+2) -= tau*tmp; 5187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp); 5197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0); 5217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Z_{k2} to annihilate T(k+1,k) 5237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k)); 5247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheRight(k+1,k,G); 5257faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.applyOnTheRight(k+1,k,G); 5267faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update Z 5277faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 5287faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(k+1,k,G.adjoint()); 5297faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(k+1,k) = Scalar(0.0); 5307faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5317faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // update x,y,z 5327faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez x = m_S.coeff(k+1,k); 5337faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez y = m_S.coeff(k+2,k); 5347faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (k < l-2) 5357faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez z = m_S.coeff(k+3,k); 5367faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // loop over k 5377faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5387faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Q_{n-1} to annihilate y = S(l,l-2) 5397faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(x,y); 5407faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheLeft(l-1,l,G.adjoint()); 5417faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.applyOnTheLeft(l-1,l,G.adjoint()); 5427faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 5437faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Q.applyOnTheRight(l-1,l,G); 5447faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.coeffRef(l,l-2) = Scalar(0.0); 5457faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5467faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // Z_{n-1} to annihilate T(l,l-1) 5477faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1)); 5487faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S.applyOnTheRight(l,l-1,G); 5497faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.applyOnTheRight(l,l-1,G); 5507faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (m_computeQZ) 5517faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_Z.applyOnTheLeft(l,l-1,G.adjoint()); 5527faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_T.coeffRef(l,l-1) = Scalar(0.0); 5537faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5547faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5557faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5567faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez template<typename MatrixType> 5577faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) 5587faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5597faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5607faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez const Index dim = A_in.cols(); 5617faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5627faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert (A_in.rows()==dim && A_in.cols()==dim 5637faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && B_in.rows()==dim && B_in.cols()==dim 5647faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez && "Need square matrices of the same dimension"); 5657faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5667faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_isInitialized = true; 5677faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_computeQZ = computeQZ; 5687faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_S = A_in; m_T = B_in; 5697faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_workspace.resize(dim*2); 5707faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_global_iter = 0; 5717faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5727faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // entrance point: hessenberg triangular decomposition 5737faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez hessenbergTriangular(); 5747faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // compute L1 vector norms of T, S into m_normOfS, m_normOfT 5757faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez computeNorms(); 5767faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5777faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index l = dim-1, 5787faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez f, 5797faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez local_iter = 0; 5807faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 5817faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez while (l>0 && local_iter<m_maxIters) 5827faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5837faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez f = findSmallSubdiagEntry(l); 5847faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // now rows and columns f..l (including) decouple from the rest of the problem 5857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0); 5867faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (f == l) // One root found 5877faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5887faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez l--; 5897faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez local_iter = 0; 5907faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5917faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else if (f == l-1) // Two roots found 5927faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5937faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez splitOffTwoRows(f); 5947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez l -= 2; 5957faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez local_iter = 0; 5967faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 5977faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else // No convergence yet 5987faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 5997faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations 6007faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez Index z = findSmallDiagEntry(f,l); 6017faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez if (z>=f) 6027faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6037faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // zero found 6047faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez pushDownZero(z,f,l); 6057faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6067faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez else 6077faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez { 6087faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg 6097faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to 6107faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // apply a QR-like iteration to rows and columns f..l. 6117faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez step(f,l, local_iter); 6127faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez local_iter++; 6137faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_global_iter++; 6147faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6157faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } 6177faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez // check if we converged before reaching iterations limit 6187faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez m_info = (local_iter<m_maxIters) ? Success : NoConvergence; 6197faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez return *this; 6207faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez } // end compute 6217faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6227faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez} // end namespace Eigen 6237faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez 6247faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez#endif //EIGEN_REAL_QZ 625