| // This file is part of Eigen, a lightweight C++ template library | |
| // for linear algebra. | |
| // | |
| // Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru> | |
| // | |
| // This Source Code Form is subject to the terms of the Mozilla | |
| // Public License v. 2.0. If a copy of the MPL was not distributed | |
| // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. | |
| #ifndef EIGEN_REAL_QZ_H | |
| #define EIGEN_REAL_QZ_H | |
| namespace Eigen { | |
| /** \eigenvalues_module \ingroup Eigenvalues_Module | |
| * | |
| * | |
| * \class RealQZ | |
| * | |
| * \brief Performs a real QZ decomposition of a pair of square matrices | |
| * | |
| * \tparam _MatrixType the type of the matrix of which we are computing the | |
| * real QZ decomposition; this is expected to be an instantiation of the | |
| * Matrix class template. | |
| * | |
| * Given a real square matrices A and B, this class computes the real QZ | |
| * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are | |
| * real orthogonal matrixes, T is upper-triangular matrix, and S is upper | |
| * quasi-triangular matrix. An orthogonal matrix is a matrix whose | |
| * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular | |
| * matrix is a block-triangular matrix whose diagonal consists of 1-by-1 | |
| * blocks and 2-by-2 blocks where further reduction is impossible due to | |
| * complex eigenvalues. | |
| * | |
| * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from | |
| * 1x1 and 2x2 blocks on the diagonals of S and T. | |
| * | |
| * Call the function compute() to compute the real QZ decomposition of a | |
| * given pair of matrices. Alternatively, you can use the | |
| * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ) | |
| * constructor which computes the real QZ decomposition at construction | |
| * time. Once the decomposition is computed, you can use the matrixS(), | |
| * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices | |
| * S, T, Q and Z in the decomposition. If computeQZ==false, some time | |
| * is saved by not computing matrices Q and Z. | |
| * | |
| * Example: \include RealQZ_compute.cpp | |
| * Output: \include RealQZ_compute.out | |
| * | |
| * \note The implementation is based on the algorithm in "Matrix Computations" | |
| * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for | |
| * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart. | |
| * | |
| * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver | |
| */ | |
| template<typename _MatrixType> class RealQZ | |
| { | |
| public: | |
| typedef _MatrixType MatrixType; | |
| enum { | |
| RowsAtCompileTime = MatrixType::RowsAtCompileTime, | |
| ColsAtCompileTime = MatrixType::ColsAtCompileTime, | |
| Options = MatrixType::Options, | |
| MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, | |
| MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime | |
| }; | |
| typedef typename MatrixType::Scalar Scalar; | |
| typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; | |
| typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 | |
| typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; | |
| typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; | |
| /** \brief Default constructor. | |
| * | |
| * \param [in] size Positive integer, size of the matrix whose QZ decomposition will be computed. | |
| * | |
| * The default constructor is useful in cases in which the user intends to | |
| * perform decompositions via compute(). The \p size parameter is only | |
| * used as a hint. It is not an error to give a wrong \p size, but it may | |
| * impair performance. | |
| * | |
| * \sa compute() for an example. | |
| */ | |
| explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : | |
| m_S(size, size), | |
| m_T(size, size), | |
| m_Q(size, size), | |
| m_Z(size, size), | |
| m_workspace(size*2), | |
| m_maxIters(400), | |
| m_isInitialized(false), | |
| m_computeQZ(true) | |
| {} | |
| /** \brief Constructor; computes real QZ decomposition of given matrices | |
| * | |
| * \param[in] A Matrix A. | |
| * \param[in] B Matrix B. | |
| * \param[in] computeQZ If false, A and Z are not computed. | |
| * | |
| * This constructor calls compute() to compute the QZ decomposition. | |
| */ | |
| RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) : | |
| m_S(A.rows(),A.cols()), | |
| m_T(A.rows(),A.cols()), | |
| m_Q(A.rows(),A.cols()), | |
| m_Z(A.rows(),A.cols()), | |
| m_workspace(A.rows()*2), | |
| m_maxIters(400), | |
| m_isInitialized(false), | |
| m_computeQZ(true) | |
| { | |
| compute(A, B, computeQZ); | |
| } | |
| /** \brief Returns matrix Q in the QZ decomposition. | |
| * | |
| * \returns A const reference to the matrix Q. | |
| */ | |
| const MatrixType& matrixQ() const { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); | |
| return m_Q; | |
| } | |
| /** \brief Returns matrix Z in the QZ decomposition. | |
| * | |
| * \returns A const reference to the matrix Z. | |
| */ | |
| const MatrixType& matrixZ() const { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); | |
| return m_Z; | |
| } | |
| /** \brief Returns matrix S in the QZ decomposition. | |
| * | |
| * \returns A const reference to the matrix S. | |
| */ | |
| const MatrixType& matrixS() const { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| return m_S; | |
| } | |
| /** \brief Returns matrix S in the QZ decomposition. | |
| * | |
| * \returns A const reference to the matrix S. | |
| */ | |
| const MatrixType& matrixT() const { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| return m_T; | |
| } | |
| /** \brief Computes QZ decomposition of given matrix. | |
| * | |
| * \param[in] A Matrix A. | |
| * \param[in] B Matrix B. | |
| * \param[in] computeQZ If false, A and Z are not computed. | |
| * \returns Reference to \c *this | |
| */ | |
| RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true); | |
| /** \brief Reports whether previous computation was successful. | |
| * | |
| * \returns \c Success if computation was successful, \c NoConvergence otherwise. | |
| */ | |
| ComputationInfo info() const | |
| { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| return m_info; | |
| } | |
| /** \brief Returns number of performed QR-like iterations. | |
| */ | |
| Index iterations() const | |
| { | |
| eigen_assert(m_isInitialized && "RealQZ is not initialized."); | |
| return m_global_iter; | |
| } | |
| /** Sets the maximal number of iterations allowed to converge to one eigenvalue | |
| * or decouple the problem. | |
| */ | |
| RealQZ& setMaxIterations(Index maxIters) | |
| { | |
| m_maxIters = maxIters; | |
| return *this; | |
| } | |
| private: | |
| MatrixType m_S, m_T, m_Q, m_Z; | |
| Matrix<Scalar,Dynamic,1> m_workspace; | |
| ComputationInfo m_info; | |
| Index m_maxIters; | |
| bool m_isInitialized; | |
| bool m_computeQZ; | |
| Scalar m_normOfT, m_normOfS; | |
| Index m_global_iter; | |
| typedef Matrix<Scalar,3,1> Vector3s; | |
| typedef Matrix<Scalar,2,1> Vector2s; | |
| typedef Matrix<Scalar,2,2> Matrix2s; | |
| typedef JacobiRotation<Scalar> JRs; | |
| void hessenbergTriangular(); | |
| void computeNorms(); | |
| Index findSmallSubdiagEntry(Index iu); | |
| Index findSmallDiagEntry(Index f, Index l); | |
| void splitOffTwoRows(Index i); | |
| void pushDownZero(Index z, Index f, Index l); | |
| void step(Index f, Index l, Index iter); | |
| }; // RealQZ | |
| /** \internal Reduces S and T to upper Hessenberg - triangular form */ | |
| template<typename MatrixType> | |
| void RealQZ<MatrixType>::hessenbergTriangular() | |
| { | |
| const Index dim = m_S.cols(); | |
| // perform QR decomposition of T, overwrite T with R, save Q | |
| HouseholderQR<MatrixType> qrT(m_T); | |
| m_T = qrT.matrixQR(); | |
| m_T.template triangularView<StrictlyLower>().setZero(); | |
| m_Q = qrT.householderQ(); | |
| // overwrite S with Q* S | |
| m_S.applyOnTheLeft(m_Q.adjoint()); | |
| // init Z as Identity | |
| if (m_computeQZ) | |
| m_Z = MatrixType::Identity(dim,dim); | |
| // reduce S to upper Hessenberg with Givens rotations | |
| for (Index j=0; j<=dim-3; j++) { | |
| for (Index i=dim-1; i>=j+2; i--) { | |
| JRs G; | |
| // kill S(i,j) | |
| if(m_S.coeff(i,j) != 0) | |
| { | |
| G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j)); | |
| m_S.coeffRef(i,j) = Scalar(0.0); | |
| m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint()); | |
| m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint()); | |
| // update Q | |
| if (m_computeQZ) | |
| m_Q.applyOnTheRight(i-1,i,G); | |
| } | |
| // kill T(i,i-1) | |
| if(m_T.coeff(i,i-1)!=Scalar(0)) | |
| { | |
| G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i)); | |
| m_T.coeffRef(i,i-1) = Scalar(0.0); | |
| m_S.applyOnTheRight(i,i-1,G); | |
| m_T.topRows(i).applyOnTheRight(i,i-1,G); | |
| // update Z | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(i,i-1,G.adjoint()); | |
| } | |
| } | |
| } | |
| } | |
| /** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */ | |
| template<typename MatrixType> | |
| inline void RealQZ<MatrixType>::computeNorms() | |
| { | |
| const Index size = m_S.cols(); | |
| m_normOfS = Scalar(0.0); | |
| m_normOfT = Scalar(0.0); | |
| for (Index j = 0; j < size; ++j) | |
| { | |
| m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); | |
| m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum(); | |
| } | |
| } | |
| /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ | |
| template<typename MatrixType> | |
| inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) | |
| { | |
| using std::abs; | |
| Index res = iu; | |
| while (res > 0) | |
| { | |
| Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res)); | |
| if (s == Scalar(0.0)) | |
| s = m_normOfS; | |
| if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s) | |
| break; | |
| res--; | |
| } | |
| return res; | |
| } | |
| /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */ | |
| template<typename MatrixType> | |
| inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) | |
| { | |
| using std::abs; | |
| Index res = l; | |
| while (res >= f) { | |
| if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT) | |
| break; | |
| res--; | |
| } | |
| return res; | |
| } | |
| /** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */ | |
| template<typename MatrixType> | |
| inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) | |
| { | |
| using std::abs; | |
| using std::sqrt; | |
| const Index dim=m_S.cols(); | |
| if (abs(m_S.coeff(i+1,i))==Scalar(0)) | |
| return; | |
| Index j = findSmallDiagEntry(i,i+1); | |
| if (j==i-1) | |
| { | |
| // block of (S T^{-1}) | |
| Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>(). | |
| template solve<OnTheRight>(m_S.template block<2,2>(i,i)); | |
| Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1)); | |
| Scalar q = p*p + STi(1,0)*STi(0,1); | |
| if (q>=0) { | |
| Scalar z = sqrt(q); | |
| // one QR-like iteration for ABi - lambda I | |
| // is enough - when we know exact eigenvalue in advance, | |
| // convergence is immediate | |
| JRs G; | |
| if (p>=0) | |
| G.makeGivens(p + z, STi(1,0)); | |
| else | |
| G.makeGivens(p - z, STi(1,0)); | |
| m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); | |
| m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); | |
| // update Q | |
| if (m_computeQZ) | |
| m_Q.applyOnTheRight(i,i+1,G); | |
| G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i)); | |
| m_S.topRows(i+2).applyOnTheRight(i+1,i,G); | |
| m_T.topRows(i+2).applyOnTheRight(i+1,i,G); | |
| // update Z | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(i+1,i,G.adjoint()); | |
| m_S.coeffRef(i+1,i) = Scalar(0.0); | |
| m_T.coeffRef(i+1,i) = Scalar(0.0); | |
| } | |
| } | |
| else | |
| { | |
| pushDownZero(j,i,i+1); | |
| } | |
| } | |
| /** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */ | |
| template<typename MatrixType> | |
| inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l) | |
| { | |
| JRs G; | |
| const Index dim = m_S.cols(); | |
| for (Index zz=z; zz<l; zz++) | |
| { | |
| // push 0 down | |
| Index firstColS = zz>f ? (zz-1) : zz; | |
| G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1)); | |
| m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint()); | |
| m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint()); | |
| m_T.coeffRef(zz+1,zz+1) = Scalar(0.0); | |
| // update Q | |
| if (m_computeQZ) | |
| m_Q.applyOnTheRight(zz,zz+1,G); | |
| // kill S(zz+1, zz-1) | |
| if (zz>f) | |
| { | |
| G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1)); | |
| m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G); | |
| m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G); | |
| m_S.coeffRef(zz+1,zz-1) = Scalar(0.0); | |
| // update Z | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(zz,zz-1,G.adjoint()); | |
| } | |
| } | |
| // finally kill S(l,l-1) | |
| G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1)); | |
| m_S.applyOnTheRight(l,l-1,G); | |
| m_T.applyOnTheRight(l,l-1,G); | |
| m_S.coeffRef(l,l-1)=Scalar(0.0); | |
| // update Z | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(l,l-1,G.adjoint()); | |
| } | |
| /** \internal QR-like iterative step for block f..l */ | |
| template<typename MatrixType> | |
| inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) | |
| { | |
| using std::abs; | |
| const Index dim = m_S.cols(); | |
| // x, y, z | |
| Scalar x, y, z; | |
| if (iter==10) | |
| { | |
| // Wilkinson ad hoc shift | |
| const Scalar | |
| a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1), | |
| a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1), | |
| b12=m_T.coeff(f+0,f+1), | |
| b11i=Scalar(1.0)/m_T.coeff(f+0,f+0), | |
| b22i=Scalar(1.0)/m_T.coeff(f+1,f+1), | |
| a87=m_S.coeff(l-1,l-2), | |
| a98=m_S.coeff(l-0,l-1), | |
| b77i=Scalar(1.0)/m_T.coeff(l-2,l-2), | |
| b88i=Scalar(1.0)/m_T.coeff(l-1,l-1); | |
| Scalar ss = abs(a87*b77i) + abs(a98*b88i), | |
| lpl = Scalar(1.5)*ss, | |
| ll = ss*ss; | |
| x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i | |
| - a11*a21*b12*b11i*b11i*b22i; | |
| y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i | |
| - a21*a21*b12*b11i*b11i*b22i; | |
| z = a21*a32*b11i*b22i; | |
| } | |
| else if (iter==16) | |
| { | |
| // another exceptional shift | |
| 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) / | |
| (m_T.coeff(l-1,l-1)*m_T.coeff(l,l)); | |
| y = m_S.coeff(f+1,f)/m_T.coeff(f,f); | |
| z = 0; | |
| } | |
| else if (iter>23 && !(iter%8)) | |
| { | |
| // extremely exceptional shift | |
| x = internal::random<Scalar>(-1.0,1.0); | |
| y = internal::random<Scalar>(-1.0,1.0); | |
| z = internal::random<Scalar>(-1.0,1.0); | |
| } | |
| else | |
| { | |
| // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1 | |
| // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where | |
| // U and V are 2x2 bottom right sub matrices of A and B. Thus: | |
| // = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1) | |
| // = AB^-1AB^-1 + det(M) - tr(M)(AB^-1) | |
| // Since we are only interested in having x, y, z with a correct ratio, we have: | |
| const Scalar | |
| a11 = m_S.coeff(f,f), a12 = m_S.coeff(f,f+1), | |
| a21 = m_S.coeff(f+1,f), a22 = m_S.coeff(f+1,f+1), | |
| a32 = m_S.coeff(f+2,f+1), | |
| a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l), | |
| a98 = m_S.coeff(l,l-1), a99 = m_S.coeff(l,l), | |
| b11 = m_T.coeff(f,f), b12 = m_T.coeff(f,f+1), | |
| b22 = m_T.coeff(f+1,f+1), | |
| b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l), | |
| b99 = m_T.coeff(l,l); | |
| x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21) | |
| + a12/b22 - (a11/b11)*(b12/b22); | |
| y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99); | |
| z = a32/b22; | |
| } | |
| JRs G; | |
| for (Index k=f; k<=l-2; k++) | |
| { | |
| // variables for Householder reflections | |
| Vector2s essential2; | |
| Scalar tau, beta; | |
| Vector3s hr(x,y,z); | |
| // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1) | |
| hr.makeHouseholderInPlace(tau, beta); | |
| essential2 = hr.template bottomRows<2>(); | |
| Index fc=(std::max)(k-1,Index(0)); // first col to update | |
| m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); | |
| m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); | |
| if (m_computeQZ) | |
| m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data()); | |
| if (k>f) | |
| m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0); | |
| // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k) | |
| hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1); | |
| hr.makeHouseholderInPlace(tau, beta); | |
| essential2 = hr.template bottomRows<2>(); | |
| { | |
| Index lr = (std::min)(k+4,dim); // last row to update | |
| Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr); | |
| // S | |
| tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2; | |
| tmp += m_S.col(k+2).head(lr); | |
| m_S.col(k+2).head(lr) -= tau*tmp; | |
| m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); | |
| // T | |
| tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2; | |
| tmp += m_T.col(k+2).head(lr); | |
| m_T.col(k+2).head(lr) -= tau*tmp; | |
| m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); | |
| } | |
| if (m_computeQZ) | |
| { | |
| // Z | |
| Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim); | |
| tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k)); | |
| tmp += m_Z.row(k+2); | |
| m_Z.row(k+2) -= tau*tmp; | |
| m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp); | |
| } | |
| m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0); | |
| // Z_{k2} to annihilate T(k+1,k) | |
| G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k)); | |
| m_S.applyOnTheRight(k+1,k,G); | |
| m_T.applyOnTheRight(k+1,k,G); | |
| // update Z | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(k+1,k,G.adjoint()); | |
| m_T.coeffRef(k+1,k) = Scalar(0.0); | |
| // update x,y,z | |
| x = m_S.coeff(k+1,k); | |
| y = m_S.coeff(k+2,k); | |
| if (k < l-2) | |
| z = m_S.coeff(k+3,k); | |
| } // loop over k | |
| // Q_{n-1} to annihilate y = S(l,l-2) | |
| G.makeGivens(x,y); | |
| m_S.applyOnTheLeft(l-1,l,G.adjoint()); | |
| m_T.applyOnTheLeft(l-1,l,G.adjoint()); | |
| if (m_computeQZ) | |
| m_Q.applyOnTheRight(l-1,l,G); | |
| m_S.coeffRef(l,l-2) = Scalar(0.0); | |
| // Z_{n-1} to annihilate T(l,l-1) | |
| G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1)); | |
| m_S.applyOnTheRight(l,l-1,G); | |
| m_T.applyOnTheRight(l,l-1,G); | |
| if (m_computeQZ) | |
| m_Z.applyOnTheLeft(l,l-1,G.adjoint()); | |
| m_T.coeffRef(l,l-1) = Scalar(0.0); | |
| } | |
| template<typename MatrixType> | |
| RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) | |
| { | |
| const Index dim = A_in.cols(); | |
| eigen_assert (A_in.rows()==dim && A_in.cols()==dim | |
| && B_in.rows()==dim && B_in.cols()==dim | |
| && "Need square matrices of the same dimension"); | |
| m_isInitialized = true; | |
| m_computeQZ = computeQZ; | |
| m_S = A_in; m_T = B_in; | |
| m_workspace.resize(dim*2); | |
| m_global_iter = 0; | |
| // entrance point: hessenberg triangular decomposition | |
| hessenbergTriangular(); | |
| // compute L1 vector norms of T, S into m_normOfS, m_normOfT | |
| computeNorms(); | |
| Index l = dim-1, | |
| f, | |
| local_iter = 0; | |
| while (l>0 && local_iter<m_maxIters) | |
| { | |
| f = findSmallSubdiagEntry(l); | |
| // now rows and columns f..l (including) decouple from the rest of the problem | |
| if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0); | |
| if (f == l) // One root found | |
| { | |
| l--; | |
| local_iter = 0; | |
| } | |
| else if (f == l-1) // Two roots found | |
| { | |
| splitOffTwoRows(f); | |
| l -= 2; | |
| local_iter = 0; | |
| } | |
| else // No convergence yet | |
| { | |
| // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations | |
| Index z = findSmallDiagEntry(f,l); | |
| if (z>=f) | |
| { | |
| // zero found | |
| pushDownZero(z,f,l); | |
| } | |
| else | |
| { | |
| // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg | |
| // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to | |
| // apply a QR-like iteration to rows and columns f..l. | |
| step(f,l, local_iter); | |
| local_iter++; | |
| m_global_iter++; | |
| } | |
| } | |
| } | |
| // check if we converged before reaching iterations limit | |
| m_info = (local_iter<m_maxIters) ? Success : NoConvergence; | |
| // For each non triangular 2x2 diagonal block of S, | |
| // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD. | |
| // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors, | |
| // and is in par with Lapack/Matlab QZ. | |
| if(m_info==Success) | |
| { | |
| for(Index i=0; i<dim-1; ++i) | |
| { | |
| if(m_S.coeff(i+1, i) != Scalar(0)) | |
| { | |
| JacobiRotation<Scalar> j_left, j_right; | |
| internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right); | |
| // Apply resulting Jacobi rotations | |
| m_S.applyOnTheLeft(i,i+1,j_left); | |
| m_S.applyOnTheRight(i,i+1,j_right); | |
| m_T.applyOnTheLeft(i,i+1,j_left); | |
| m_T.applyOnTheRight(i,i+1,j_right); | |
| m_T(i+1,i) = m_T(i,i+1) = Scalar(0); | |
| if(m_computeQZ) { | |
| m_Q.applyOnTheRight(i,i+1,j_left.transpose()); | |
| m_Z.applyOnTheLeft(i,i+1,j_right.transpose()); | |
| } | |
| i++; | |
| } | |
| } | |
| } | |
| return *this; | |
| } // end compute | |
| } // end namespace Eigen | |
| #endif //EIGEN_REAL_QZ |