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