1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_REAL_QZ_H 11 #define EIGEN_REAL_QZ_H 12 13 namespace Eigen { 14 15 /** \eigenvalues_module \ingroup Eigenvalues_Module 16 * 17 * 18 * \class RealQZ 19 * 20 * \brief Performs a real QZ decomposition of a pair of square matrices 21 * 22 * \tparam _MatrixType the type of the matrix of which we are computing the 23 * real QZ decomposition; this is expected to be an instantiation of the 24 * Matrix class template. 25 * 26 * Given a real square matrices A and B, this class computes the real QZ 27 * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are 28 * real orthogonal matrixes, T is upper-triangular matrix, and S is upper 29 * quasi-triangular matrix. An orthogonal matrix is a matrix whose 30 * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular 31 * matrix is a block-triangular matrix whose diagonal consists of 1-by-1 32 * blocks and 2-by-2 blocks where further reduction is impossible due to 33 * complex eigenvalues. 34 * 35 * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from 36 * 1x1 and 2x2 blocks on the diagonals of S and T. 37 * 38 * Call the function compute() to compute the real QZ decomposition of a 39 * given pair of matrices. Alternatively, you can use the 40 * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ) 41 * constructor which computes the real QZ decomposition at construction 42 * time. Once the decomposition is computed, you can use the matrixS(), 43 * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices 44 * S, T, Q and Z in the decomposition. If computeQZ==false, some time 45 * is saved by not computing matrices Q and Z. 46 * 47 * Example: \include RealQZ_compute.cpp 48 * Output: \include RealQZ_compute.out 49 * 50 * \note The implementation is based on the algorithm in "Matrix Computations" 51 * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for 52 * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart. 53 * 54 * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver 55 */ 56 57 template<typename _MatrixType> class RealQZ 58 { 59 public: 60 typedef _MatrixType MatrixType; 61 enum { 62 RowsAtCompileTime = MatrixType::RowsAtCompileTime, 63 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 64 Options = MatrixType::Options, 65 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 66 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 67 }; 68 typedef typename MatrixType::Scalar Scalar; 69 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; 70 typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3 71 72 typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType; 73 typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType; 74 75 /** \brief Default constructor. 76 * 77 * \param [in] size Positive integer, size of the matrix whose QZ decomposition will be computed. 78 * 79 * The default constructor is useful in cases in which the user intends to 80 * perform decompositions via compute(). The \p size parameter is only 81 * used as a hint. It is not an error to give a wrong \p size, but it may 82 * impair performance. 83 * 84 * \sa compute() for an example. 85 */ 86 explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_S(size,size)87 m_S(size, size), 88 m_T(size, size), 89 m_Q(size, size), 90 m_Z(size, size), 91 m_workspace(size*2), 92 m_maxIters(400), 93 m_isInitialized(false) 94 { } 95 96 /** \brief Constructor; computes real QZ decomposition of given matrices 97 * 98 * \param[in] A Matrix A. 99 * \param[in] B Matrix B. 100 * \param[in] computeQZ If false, A and Z are not computed. 101 * 102 * This constructor calls compute() to compute the QZ decomposition. 103 */ 104 RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) : 105 m_S(A.rows(),A.cols()), 106 m_T(A.rows(),A.cols()), 107 m_Q(A.rows(),A.cols()), 108 m_Z(A.rows(),A.cols()), 109 m_workspace(A.rows()*2), 110 m_maxIters(400), 111 m_isInitialized(false) { 112 compute(A, B, computeQZ); 113 } 114 115 /** \brief Returns matrix Q in the QZ decomposition. 116 * 117 * \returns A const reference to the matrix Q. 118 */ matrixQ()119 const MatrixType& matrixQ() const { 120 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 121 eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); 122 return m_Q; 123 } 124 125 /** \brief Returns matrix Z in the QZ decomposition. 126 * 127 * \returns A const reference to the matrix Z. 128 */ matrixZ()129 const MatrixType& matrixZ() const { 130 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 131 eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition."); 132 return m_Z; 133 } 134 135 /** \brief Returns matrix S in the QZ decomposition. 136 * 137 * \returns A const reference to the matrix S. 138 */ matrixS()139 const MatrixType& matrixS() const { 140 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 141 return m_S; 142 } 143 144 /** \brief Returns matrix S in the QZ decomposition. 145 * 146 * \returns A const reference to the matrix S. 147 */ matrixT()148 const MatrixType& matrixT() const { 149 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 150 return m_T; 151 } 152 153 /** \brief Computes QZ decomposition of given matrix. 154 * 155 * \param[in] A Matrix A. 156 * \param[in] B Matrix B. 157 * \param[in] computeQZ If false, A and Z are not computed. 158 * \returns Reference to \c *this 159 */ 160 RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true); 161 162 /** \brief Reports whether previous computation was successful. 163 * 164 * \returns \c Success if computation was succesful, \c NoConvergence otherwise. 165 */ info()166 ComputationInfo info() const 167 { 168 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 169 return m_info; 170 } 171 172 /** \brief Returns number of performed QR-like iterations. 173 */ iterations()174 Index iterations() const 175 { 176 eigen_assert(m_isInitialized && "RealQZ is not initialized."); 177 return m_global_iter; 178 } 179 180 /** Sets the maximal number of iterations allowed to converge to one eigenvalue 181 * or decouple the problem. 182 */ setMaxIterations(Index maxIters)183 RealQZ& setMaxIterations(Index maxIters) 184 { 185 m_maxIters = maxIters; 186 return *this; 187 } 188 189 private: 190 191 MatrixType m_S, m_T, m_Q, m_Z; 192 Matrix<Scalar,Dynamic,1> m_workspace; 193 ComputationInfo m_info; 194 Index m_maxIters; 195 bool m_isInitialized; 196 bool m_computeQZ; 197 Scalar m_normOfT, m_normOfS; 198 Index m_global_iter; 199 200 typedef Matrix<Scalar,3,1> Vector3s; 201 typedef Matrix<Scalar,2,1> Vector2s; 202 typedef Matrix<Scalar,2,2> Matrix2s; 203 typedef JacobiRotation<Scalar> JRs; 204 205 void hessenbergTriangular(); 206 void computeNorms(); 207 Index findSmallSubdiagEntry(Index iu); 208 Index findSmallDiagEntry(Index f, Index l); 209 void splitOffTwoRows(Index i); 210 void pushDownZero(Index z, Index f, Index l); 211 void step(Index f, Index l, Index iter); 212 213 }; // RealQZ 214 215 /** \internal Reduces S and T to upper Hessenberg - triangular form */ 216 template<typename MatrixType> hessenbergTriangular()217 void RealQZ<MatrixType>::hessenbergTriangular() 218 { 219 220 const Index dim = m_S.cols(); 221 222 // perform QR decomposition of T, overwrite T with R, save Q 223 HouseholderQR<MatrixType> qrT(m_T); 224 m_T = qrT.matrixQR(); 225 m_T.template triangularView<StrictlyLower>().setZero(); 226 m_Q = qrT.householderQ(); 227 // overwrite S with Q* S 228 m_S.applyOnTheLeft(m_Q.adjoint()); 229 // init Z as Identity 230 if (m_computeQZ) 231 m_Z = MatrixType::Identity(dim,dim); 232 // reduce S to upper Hessenberg with Givens rotations 233 for (Index j=0; j<=dim-3; j++) { 234 for (Index i=dim-1; i>=j+2; i--) { 235 JRs G; 236 // kill S(i,j) 237 if(m_S.coeff(i,j) != 0) 238 { 239 G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j)); 240 m_S.coeffRef(i,j) = Scalar(0.0); 241 m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint()); 242 m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint()); 243 // update Q 244 if (m_computeQZ) 245 m_Q.applyOnTheRight(i-1,i,G); 246 } 247 // kill T(i,i-1) 248 if(m_T.coeff(i,i-1)!=Scalar(0)) 249 { 250 G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i)); 251 m_T.coeffRef(i,i-1) = Scalar(0.0); 252 m_S.applyOnTheRight(i,i-1,G); 253 m_T.topRows(i).applyOnTheRight(i,i-1,G); 254 // update Z 255 if (m_computeQZ) 256 m_Z.applyOnTheLeft(i,i-1,G.adjoint()); 257 } 258 } 259 } 260 } 261 262 /** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */ 263 template<typename MatrixType> computeNorms()264 inline void RealQZ<MatrixType>::computeNorms() 265 { 266 const Index size = m_S.cols(); 267 m_normOfS = Scalar(0.0); 268 m_normOfT = Scalar(0.0); 269 for (Index j = 0; j < size; ++j) 270 { 271 m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum(); 272 m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum(); 273 } 274 } 275 276 277 /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */ 278 template<typename MatrixType> findSmallSubdiagEntry(Index iu)279 inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) 280 { 281 using std::abs; 282 Index res = iu; 283 while (res > 0) 284 { 285 Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res)); 286 if (s == Scalar(0.0)) 287 s = m_normOfS; 288 if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s) 289 break; 290 res--; 291 } 292 return res; 293 } 294 295 /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1) */ 296 template<typename MatrixType> findSmallDiagEntry(Index f,Index l)297 inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) 298 { 299 using std::abs; 300 Index res = l; 301 while (res >= f) { 302 if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT) 303 break; 304 res--; 305 } 306 return res; 307 } 308 309 /** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */ 310 template<typename MatrixType> splitOffTwoRows(Index i)311 inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) 312 { 313 using std::abs; 314 using std::sqrt; 315 const Index dim=m_S.cols(); 316 if (abs(m_S.coeff(i+1,i))==Scalar(0)) 317 return; 318 Index j = findSmallDiagEntry(i,i+1); 319 if (j==i-1) 320 { 321 // block of (S T^{-1}) 322 Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>(). 323 template solve<OnTheRight>(m_S.template block<2,2>(i,i)); 324 Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1)); 325 Scalar q = p*p + STi(1,0)*STi(0,1); 326 if (q>=0) { 327 Scalar z = sqrt(q); 328 // one QR-like iteration for ABi - lambda I 329 // is enough - when we know exact eigenvalue in advance, 330 // convergence is immediate 331 JRs G; 332 if (p>=0) 333 G.makeGivens(p + z, STi(1,0)); 334 else 335 G.makeGivens(p - z, STi(1,0)); 336 m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); 337 m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint()); 338 // update Q 339 if (m_computeQZ) 340 m_Q.applyOnTheRight(i,i+1,G); 341 342 G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i)); 343 m_S.topRows(i+2).applyOnTheRight(i+1,i,G); 344 m_T.topRows(i+2).applyOnTheRight(i+1,i,G); 345 // update Z 346 if (m_computeQZ) 347 m_Z.applyOnTheLeft(i+1,i,G.adjoint()); 348 349 m_S.coeffRef(i+1,i) = Scalar(0.0); 350 m_T.coeffRef(i+1,i) = Scalar(0.0); 351 } 352 } 353 else 354 { 355 pushDownZero(j,i,i+1); 356 } 357 } 358 359 /** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */ 360 template<typename MatrixType> pushDownZero(Index z,Index f,Index l)361 inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l) 362 { 363 JRs G; 364 const Index dim = m_S.cols(); 365 for (Index zz=z; zz<l; zz++) 366 { 367 // push 0 down 368 Index firstColS = zz>f ? (zz-1) : zz; 369 G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1)); 370 m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint()); 371 m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint()); 372 m_T.coeffRef(zz+1,zz+1) = Scalar(0.0); 373 // update Q 374 if (m_computeQZ) 375 m_Q.applyOnTheRight(zz,zz+1,G); 376 // kill S(zz+1, zz-1) 377 if (zz>f) 378 { 379 G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1)); 380 m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G); 381 m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G); 382 m_S.coeffRef(zz+1,zz-1) = Scalar(0.0); 383 // update Z 384 if (m_computeQZ) 385 m_Z.applyOnTheLeft(zz,zz-1,G.adjoint()); 386 } 387 } 388 // finally kill S(l,l-1) 389 G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1)); 390 m_S.applyOnTheRight(l,l-1,G); 391 m_T.applyOnTheRight(l,l-1,G); 392 m_S.coeffRef(l,l-1)=Scalar(0.0); 393 // update Z 394 if (m_computeQZ) 395 m_Z.applyOnTheLeft(l,l-1,G.adjoint()); 396 } 397 398 /** \internal QR-like iterative step for block f..l */ 399 template<typename MatrixType> step(Index f,Index l,Index iter)400 inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) 401 { 402 using std::abs; 403 const Index dim = m_S.cols(); 404 405 // x, y, z 406 Scalar x, y, z; 407 if (iter==10) 408 { 409 // Wilkinson ad hoc shift 410 const Scalar 411 a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1), 412 a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1), 413 b12=m_T.coeff(f+0,f+1), 414 b11i=Scalar(1.0)/m_T.coeff(f+0,f+0), 415 b22i=Scalar(1.0)/m_T.coeff(f+1,f+1), 416 a87=m_S.coeff(l-1,l-2), 417 a98=m_S.coeff(l-0,l-1), 418 b77i=Scalar(1.0)/m_T.coeff(l-2,l-2), 419 b88i=Scalar(1.0)/m_T.coeff(l-1,l-1); 420 Scalar ss = abs(a87*b77i) + abs(a98*b88i), 421 lpl = Scalar(1.5)*ss, 422 ll = ss*ss; 423 x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i 424 - a11*a21*b12*b11i*b11i*b22i; 425 y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i 426 - a21*a21*b12*b11i*b11i*b22i; 427 z = a21*a32*b11i*b22i; 428 } 429 else if (iter==16) 430 { 431 // another exceptional shift 432 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) / 433 (m_T.coeff(l-1,l-1)*m_T.coeff(l,l)); 434 y = m_S.coeff(f+1,f)/m_T.coeff(f,f); 435 z = 0; 436 } 437 else if (iter>23 && !(iter%8)) 438 { 439 // extremely exceptional shift 440 x = internal::random<Scalar>(-1.0,1.0); 441 y = internal::random<Scalar>(-1.0,1.0); 442 z = internal::random<Scalar>(-1.0,1.0); 443 } 444 else 445 { 446 // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1 447 // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where 448 // U and V are 2x2 bottom right sub matrices of A and B. Thus: 449 // = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1) 450 // = AB^-1AB^-1 + det(M) - tr(M)(AB^-1) 451 // Since we are only interested in having x, y, z with a correct ratio, we have: 452 const Scalar 453 a11 = m_S.coeff(f,f), a12 = m_S.coeff(f,f+1), 454 a21 = m_S.coeff(f+1,f), a22 = m_S.coeff(f+1,f+1), 455 a32 = m_S.coeff(f+2,f+1), 456 457 a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l), 458 a98 = m_S.coeff(l,l-1), a99 = m_S.coeff(l,l), 459 460 b11 = m_T.coeff(f,f), b12 = m_T.coeff(f,f+1), 461 b22 = m_T.coeff(f+1,f+1), 462 463 b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l), 464 b99 = m_T.coeff(l,l); 465 466 x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21) 467 + a12/b22 - (a11/b11)*(b12/b22); 468 y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99); 469 z = a32/b22; 470 } 471 472 JRs G; 473 474 for (Index k=f; k<=l-2; k++) 475 { 476 // variables for Householder reflections 477 Vector2s essential2; 478 Scalar tau, beta; 479 480 Vector3s hr(x,y,z); 481 482 // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1) 483 hr.makeHouseholderInPlace(tau, beta); 484 essential2 = hr.template bottomRows<2>(); 485 Index fc=(std::max)(k-1,Index(0)); // first col to update 486 m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); 487 m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data()); 488 if (m_computeQZ) 489 m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data()); 490 if (k>f) 491 m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0); 492 493 // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k) 494 hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1); 495 hr.makeHouseholderInPlace(tau, beta); 496 essential2 = hr.template bottomRows<2>(); 497 { 498 Index lr = (std::min)(k+4,dim); // last row to update 499 Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr); 500 // S 501 tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2; 502 tmp += m_S.col(k+2).head(lr); 503 m_S.col(k+2).head(lr) -= tau*tmp; 504 m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); 505 // T 506 tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2; 507 tmp += m_T.col(k+2).head(lr); 508 m_T.col(k+2).head(lr) -= tau*tmp; 509 m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint(); 510 } 511 if (m_computeQZ) 512 { 513 // Z 514 Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim); 515 tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k)); 516 tmp += m_Z.row(k+2); 517 m_Z.row(k+2) -= tau*tmp; 518 m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp); 519 } 520 m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0); 521 522 // Z_{k2} to annihilate T(k+1,k) 523 G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k)); 524 m_S.applyOnTheRight(k+1,k,G); 525 m_T.applyOnTheRight(k+1,k,G); 526 // update Z 527 if (m_computeQZ) 528 m_Z.applyOnTheLeft(k+1,k,G.adjoint()); 529 m_T.coeffRef(k+1,k) = Scalar(0.0); 530 531 // update x,y,z 532 x = m_S.coeff(k+1,k); 533 y = m_S.coeff(k+2,k); 534 if (k < l-2) 535 z = m_S.coeff(k+3,k); 536 } // loop over k 537 538 // Q_{n-1} to annihilate y = S(l,l-2) 539 G.makeGivens(x,y); 540 m_S.applyOnTheLeft(l-1,l,G.adjoint()); 541 m_T.applyOnTheLeft(l-1,l,G.adjoint()); 542 if (m_computeQZ) 543 m_Q.applyOnTheRight(l-1,l,G); 544 m_S.coeffRef(l,l-2) = Scalar(0.0); 545 546 // Z_{n-1} to annihilate T(l,l-1) 547 G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1)); 548 m_S.applyOnTheRight(l,l-1,G); 549 m_T.applyOnTheRight(l,l-1,G); 550 if (m_computeQZ) 551 m_Z.applyOnTheLeft(l,l-1,G.adjoint()); 552 m_T.coeffRef(l,l-1) = Scalar(0.0); 553 } 554 555 template<typename MatrixType> compute(const MatrixType & A_in,const MatrixType & B_in,bool computeQZ)556 RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) 557 { 558 559 const Index dim = A_in.cols(); 560 561 eigen_assert (A_in.rows()==dim && A_in.cols()==dim 562 && B_in.rows()==dim && B_in.cols()==dim 563 && "Need square matrices of the same dimension"); 564 565 m_isInitialized = true; 566 m_computeQZ = computeQZ; 567 m_S = A_in; m_T = B_in; 568 m_workspace.resize(dim*2); 569 m_global_iter = 0; 570 571 // entrance point: hessenberg triangular decomposition 572 hessenbergTriangular(); 573 // compute L1 vector norms of T, S into m_normOfS, m_normOfT 574 computeNorms(); 575 576 Index l = dim-1, 577 f, 578 local_iter = 0; 579 580 while (l>0 && local_iter<m_maxIters) 581 { 582 f = findSmallSubdiagEntry(l); 583 // now rows and columns f..l (including) decouple from the rest of the problem 584 if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0); 585 if (f == l) // One root found 586 { 587 l--; 588 local_iter = 0; 589 } 590 else if (f == l-1) // Two roots found 591 { 592 splitOffTwoRows(f); 593 l -= 2; 594 local_iter = 0; 595 } 596 else // No convergence yet 597 { 598 // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations 599 Index z = findSmallDiagEntry(f,l); 600 if (z>=f) 601 { 602 // zero found 603 pushDownZero(z,f,l); 604 } 605 else 606 { 607 // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg 608 // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to 609 // apply a QR-like iteration to rows and columns f..l. 610 step(f,l, local_iter); 611 local_iter++; 612 m_global_iter++; 613 } 614 } 615 } 616 // check if we converged before reaching iterations limit 617 m_info = (local_iter<m_maxIters) ? Success : NoConvergence; 618 619 // For each non triangular 2x2 diagonal block of S, 620 // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD. 621 // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors, 622 // and is in par with Lapack/Matlab QZ. 623 if(m_info==Success) 624 { 625 for(Index i=0; i<dim-1; ++i) 626 { 627 if(m_S.coeff(i+1, i) != Scalar(0)) 628 { 629 JacobiRotation<Scalar> j_left, j_right; 630 internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right); 631 632 // Apply resulting Jacobi rotations 633 m_S.applyOnTheLeft(i,i+1,j_left); 634 m_S.applyOnTheRight(i,i+1,j_right); 635 m_T.applyOnTheLeft(i,i+1,j_left); 636 m_T.applyOnTheRight(i,i+1,j_right); 637 m_T(i+1,i) = m_T(i,i+1) = Scalar(0); 638 639 if(m_computeQZ) { 640 m_Q.applyOnTheRight(i,i+1,j_left.transpose()); 641 m_Z.applyOnTheLeft(i,i+1,j_right.transpose()); 642 } 643 644 i++; 645 } 646 } 647 } 648 649 return *this; 650 } // end compute 651 652 } // end namespace Eigen 653 654 #endif //EIGEN_REAL_QZ 655