Back to home page

EIC code displayed by LXR
 
 

    

Warning, file /include/eigen3/Eigen/src/Eigenvalues/RealQZ.h was not indexed or was modified since last indexation (in which case cross-reference links may be missing, inaccurate or erroneous).

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009 
0010 #ifndef EIGEN_REAL_QZ_H
0011 #define EIGEN_REAL_QZ_H
0012 
0013 namespace Eigen {
0014 
0015   /** \eigenvalues_module \ingroup Eigenvalues_Module
0016    *
0017    *
0018    * \class RealQZ
0019    *
0020    * \brief Performs a real QZ decomposition of a pair of square matrices
0021    *
0022    * \tparam _MatrixType the type of the matrix of which we are computing the
0023    * real QZ decomposition; this is expected to be an instantiation of the
0024    * Matrix class template.
0025    *
0026    * Given a real square matrices A and B, this class computes the real QZ
0027    * decomposition: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ where Q and Z are
0028    * real orthogonal matrixes, T is upper-triangular matrix, and S is upper
0029    * quasi-triangular matrix. An orthogonal matrix is a matrix whose
0030    * inverse is equal to its transpose, \f$ U^{-1} = U^T \f$. A quasi-triangular
0031    * matrix is a block-triangular matrix whose diagonal consists of 1-by-1
0032    * blocks and 2-by-2 blocks where further reduction is impossible due to
0033    * complex eigenvalues. 
0034    *
0035    * The eigenvalues of the pencil \f$ A - z B \f$ can be obtained from
0036    * 1x1 and 2x2 blocks on the diagonals of S and T.
0037    *
0038    * Call the function compute() to compute the real QZ decomposition of a
0039    * given pair of matrices. Alternatively, you can use the 
0040    * RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ)
0041    * constructor which computes the real QZ decomposition at construction
0042    * time. Once the decomposition is computed, you can use the matrixS(),
0043    * matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices
0044    * S, T, Q and Z in the decomposition. If computeQZ==false, some time
0045    * is saved by not computing matrices Q and Z.
0046    *
0047    * Example: \include RealQZ_compute.cpp
0048    * Output: \include RealQZ_compute.out
0049    *
0050    * \note The implementation is based on the algorithm in "Matrix Computations"
0051    * by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for
0052    * generalized eigenvalue problems" by C.B.Moler and G.W.Stewart.
0053    *
0054    * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver
0055    */
0056 
0057   template<typename _MatrixType> class RealQZ
0058   {
0059     public:
0060       typedef _MatrixType MatrixType;
0061       enum {
0062         RowsAtCompileTime = MatrixType::RowsAtCompileTime,
0063         ColsAtCompileTime = MatrixType::ColsAtCompileTime,
0064         Options = MatrixType::Options,
0065         MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0066         MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0067       };
0068       typedef typename MatrixType::Scalar Scalar;
0069       typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
0070       typedef Eigen::Index Index; ///< \deprecated since Eigen 3.3
0071 
0072       typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
0073       typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
0074 
0075       /** \brief Default constructor.
0076        *
0077        * \param [in] size  Positive integer, size of the matrix whose QZ decomposition will be computed.
0078        *
0079        * The default constructor is useful in cases in which the user intends to
0080        * perform decompositions via compute().  The \p size parameter is only
0081        * used as a hint. It is not an error to give a wrong \p size, but it may
0082        * impair performance.
0083        *
0084        * \sa compute() for an example.
0085        */
0086       explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) :
0087         m_S(size, size),
0088         m_T(size, size),
0089         m_Q(size, size),
0090         m_Z(size, size),
0091         m_workspace(size*2),
0092         m_maxIters(400),
0093         m_isInitialized(false),
0094         m_computeQZ(true)
0095       {}
0096 
0097       /** \brief Constructor; computes real QZ decomposition of given matrices
0098        * 
0099        * \param[in]  A          Matrix A.
0100        * \param[in]  B          Matrix B.
0101        * \param[in]  computeQZ  If false, A and Z are not computed.
0102        *
0103        * This constructor calls compute() to compute the QZ decomposition.
0104        */
0105       RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) :
0106         m_S(A.rows(),A.cols()),
0107         m_T(A.rows(),A.cols()),
0108         m_Q(A.rows(),A.cols()),
0109         m_Z(A.rows(),A.cols()),
0110         m_workspace(A.rows()*2),
0111         m_maxIters(400),
0112         m_isInitialized(false),
0113         m_computeQZ(true)
0114       {
0115         compute(A, B, computeQZ);
0116       }
0117 
0118       /** \brief Returns matrix Q in the QZ decomposition. 
0119        *
0120        * \returns A const reference to the matrix Q.
0121        */
0122       const MatrixType& matrixQ() const {
0123         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0124         eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
0125         return m_Q;
0126       }
0127 
0128       /** \brief Returns matrix Z in the QZ decomposition. 
0129        *
0130        * \returns A const reference to the matrix Z.
0131        */
0132       const MatrixType& matrixZ() const {
0133         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0134         eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
0135         return m_Z;
0136       }
0137 
0138       /** \brief Returns matrix S in the QZ decomposition. 
0139        *
0140        * \returns A const reference to the matrix S.
0141        */
0142       const MatrixType& matrixS() const {
0143         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0144         return m_S;
0145       }
0146 
0147       /** \brief Returns matrix S in the QZ decomposition. 
0148        *
0149        * \returns A const reference to the matrix S.
0150        */
0151       const MatrixType& matrixT() const {
0152         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0153         return m_T;
0154       }
0155 
0156       /** \brief Computes QZ decomposition of given matrix. 
0157        * 
0158        * \param[in]  A          Matrix A.
0159        * \param[in]  B          Matrix B.
0160        * \param[in]  computeQZ  If false, A and Z are not computed.
0161        * \returns    Reference to \c *this
0162        */
0163       RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);
0164 
0165       /** \brief Reports whether previous computation was successful.
0166        *
0167        * \returns \c Success if computation was successful, \c NoConvergence otherwise.
0168        */
0169       ComputationInfo info() const
0170       {
0171         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0172         return m_info;
0173       }
0174 
0175       /** \brief Returns number of performed QR-like iterations.
0176       */
0177       Index iterations() const
0178       {
0179         eigen_assert(m_isInitialized && "RealQZ is not initialized.");
0180         return m_global_iter;
0181       }
0182 
0183       /** Sets the maximal number of iterations allowed to converge to one eigenvalue
0184        * or decouple the problem.
0185       */
0186       RealQZ& setMaxIterations(Index maxIters)
0187       {
0188         m_maxIters = maxIters;
0189         return *this;
0190       }
0191 
0192     private:
0193 
0194       MatrixType m_S, m_T, m_Q, m_Z;
0195       Matrix<Scalar,Dynamic,1> m_workspace;
0196       ComputationInfo m_info;
0197       Index m_maxIters;
0198       bool m_isInitialized;
0199       bool m_computeQZ;
0200       Scalar m_normOfT, m_normOfS;
0201       Index m_global_iter;
0202 
0203       typedef Matrix<Scalar,3,1> Vector3s;
0204       typedef Matrix<Scalar,2,1> Vector2s;
0205       typedef Matrix<Scalar,2,2> Matrix2s;
0206       typedef JacobiRotation<Scalar> JRs;
0207 
0208       void hessenbergTriangular();
0209       void computeNorms();
0210       Index findSmallSubdiagEntry(Index iu);
0211       Index findSmallDiagEntry(Index f, Index l);
0212       void splitOffTwoRows(Index i);
0213       void pushDownZero(Index z, Index f, Index l);
0214       void step(Index f, Index l, Index iter);
0215 
0216   }; // RealQZ
0217 
0218   /** \internal Reduces S and T to upper Hessenberg - triangular form */
0219   template<typename MatrixType>
0220     void RealQZ<MatrixType>::hessenbergTriangular()
0221     {
0222 
0223       const Index dim = m_S.cols();
0224 
0225       // perform QR decomposition of T, overwrite T with R, save Q
0226       HouseholderQR<MatrixType> qrT(m_T);
0227       m_T = qrT.matrixQR();
0228       m_T.template triangularView<StrictlyLower>().setZero();
0229       m_Q = qrT.householderQ();
0230       // overwrite S with Q* S
0231       m_S.applyOnTheLeft(m_Q.adjoint());
0232       // init Z as Identity
0233       if (m_computeQZ)
0234         m_Z = MatrixType::Identity(dim,dim);
0235       // reduce S to upper Hessenberg with Givens rotations
0236       for (Index j=0; j<=dim-3; j++) {
0237         for (Index i=dim-1; i>=j+2; i--) {
0238           JRs G;
0239           // kill S(i,j)
0240           if(m_S.coeff(i,j) != 0)
0241           {
0242             G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
0243             m_S.coeffRef(i,j) = Scalar(0.0);
0244             m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint());
0245             m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint());
0246             // update Q
0247             if (m_computeQZ)
0248               m_Q.applyOnTheRight(i-1,i,G);
0249           }
0250           // kill T(i,i-1)
0251           if(m_T.coeff(i,i-1)!=Scalar(0))
0252           {
0253             G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
0254             m_T.coeffRef(i,i-1) = Scalar(0.0);
0255             m_S.applyOnTheRight(i,i-1,G);
0256             m_T.topRows(i).applyOnTheRight(i,i-1,G);
0257             // update Z
0258             if (m_computeQZ)
0259               m_Z.applyOnTheLeft(i,i-1,G.adjoint());
0260           }
0261         }
0262       }
0263     }
0264 
0265   /** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
0266   template<typename MatrixType>
0267     inline void RealQZ<MatrixType>::computeNorms()
0268     {
0269       const Index size = m_S.cols();
0270       m_normOfS = Scalar(0.0);
0271       m_normOfT = Scalar(0.0);
0272       for (Index j = 0; j < size; ++j)
0273       {
0274         m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
0275         m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
0276       }
0277     }
0278 
0279 
0280   /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
0281   template<typename MatrixType>
0282     inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu)
0283     {
0284       using std::abs;
0285       Index res = iu;
0286       while (res > 0)
0287       {
0288         Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
0289         if (s == Scalar(0.0))
0290           s = m_normOfS;
0291         if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
0292           break;
0293         res--;
0294       }
0295       return res;
0296     }
0297 
0298   /** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1)  */
0299   template<typename MatrixType>
0300     inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l)
0301     {
0302       using std::abs;
0303       Index res = l;
0304       while (res >= f) {
0305         if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
0306           break;
0307         res--;
0308       }
0309       return res;
0310     }
0311 
0312   /** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */
0313   template<typename MatrixType>
0314     inline void RealQZ<MatrixType>::splitOffTwoRows(Index i)
0315     {
0316       using std::abs;
0317       using std::sqrt;
0318       const Index dim=m_S.cols();
0319       if (abs(m_S.coeff(i+1,i))==Scalar(0))
0320         return;
0321       Index j = findSmallDiagEntry(i,i+1);
0322       if (j==i-1)
0323       {
0324         // block of (S T^{-1})
0325         Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>().
0326           template solve<OnTheRight>(m_S.template block<2,2>(i,i));
0327         Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1));
0328         Scalar q = p*p + STi(1,0)*STi(0,1);
0329         if (q>=0) {
0330           Scalar z = sqrt(q);
0331           // one QR-like iteration for ABi - lambda I
0332           // is enough - when we know exact eigenvalue in advance,
0333           // convergence is immediate
0334           JRs G;
0335           if (p>=0)
0336             G.makeGivens(p + z, STi(1,0));
0337           else
0338             G.makeGivens(p - z, STi(1,0));
0339           m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
0340           m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
0341           // update Q
0342           if (m_computeQZ)
0343             m_Q.applyOnTheRight(i,i+1,G);
0344 
0345           G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i));
0346           m_S.topRows(i+2).applyOnTheRight(i+1,i,G);
0347           m_T.topRows(i+2).applyOnTheRight(i+1,i,G);
0348           // update Z
0349           if (m_computeQZ)
0350             m_Z.applyOnTheLeft(i+1,i,G.adjoint());
0351 
0352           m_S.coeffRef(i+1,i) = Scalar(0.0);
0353           m_T.coeffRef(i+1,i) = Scalar(0.0);
0354         }
0355       }
0356       else
0357       {
0358         pushDownZero(j,i,i+1);
0359       }
0360     }
0361 
0362   /** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */
0363   template<typename MatrixType>
0364     inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l)
0365     {
0366       JRs G;
0367       const Index dim = m_S.cols();
0368       for (Index zz=z; zz<l; zz++)
0369       {
0370         // push 0 down
0371         Index firstColS = zz>f ? (zz-1) : zz;
0372         G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1));
0373         m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint());
0374         m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint());
0375         m_T.coeffRef(zz+1,zz+1) = Scalar(0.0);
0376         // update Q
0377         if (m_computeQZ)
0378           m_Q.applyOnTheRight(zz,zz+1,G);
0379         // kill S(zz+1, zz-1)
0380         if (zz>f)
0381         {
0382           G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1));
0383           m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G);
0384           m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G);
0385           m_S.coeffRef(zz+1,zz-1) = Scalar(0.0);
0386           // update Z
0387           if (m_computeQZ)
0388             m_Z.applyOnTheLeft(zz,zz-1,G.adjoint());
0389         }
0390       }
0391       // finally kill S(l,l-1)
0392       G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1));
0393       m_S.applyOnTheRight(l,l-1,G);
0394       m_T.applyOnTheRight(l,l-1,G);
0395       m_S.coeffRef(l,l-1)=Scalar(0.0);
0396       // update Z
0397       if (m_computeQZ)
0398         m_Z.applyOnTheLeft(l,l-1,G.adjoint());
0399     }
0400 
0401   /** \internal QR-like iterative step for block f..l */
0402   template<typename MatrixType>
0403     inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter)
0404     {
0405       using std::abs;
0406       const Index dim = m_S.cols();
0407 
0408       // x, y, z
0409       Scalar x, y, z;
0410       if (iter==10)
0411       {
0412         // Wilkinson ad hoc shift
0413         const Scalar
0414           a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1),
0415           a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1),
0416           b12=m_T.coeff(f+0,f+1),
0417           b11i=Scalar(1.0)/m_T.coeff(f+0,f+0),
0418           b22i=Scalar(1.0)/m_T.coeff(f+1,f+1),
0419           a87=m_S.coeff(l-1,l-2),
0420           a98=m_S.coeff(l-0,l-1),
0421           b77i=Scalar(1.0)/m_T.coeff(l-2,l-2),
0422           b88i=Scalar(1.0)/m_T.coeff(l-1,l-1);
0423         Scalar ss = abs(a87*b77i) + abs(a98*b88i),
0424                lpl = Scalar(1.5)*ss,
0425                ll = ss*ss;
0426         x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i
0427           - a11*a21*b12*b11i*b11i*b22i;
0428         y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i 
0429           - a21*a21*b12*b11i*b11i*b22i;
0430         z = a21*a32*b11i*b22i;
0431       }
0432       else if (iter==16)
0433       {
0434         // another exceptional shift
0435         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) /
0436           (m_T.coeff(l-1,l-1)*m_T.coeff(l,l));
0437         y = m_S.coeff(f+1,f)/m_T.coeff(f,f);
0438         z = 0;
0439       }
0440       else if (iter>23 && !(iter%8))
0441       {
0442         // extremely exceptional shift
0443         x = internal::random<Scalar>(-1.0,1.0);
0444         y = internal::random<Scalar>(-1.0,1.0);
0445         z = internal::random<Scalar>(-1.0,1.0);
0446       }
0447       else
0448       {
0449         // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1
0450         // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where
0451         // U and V are 2x2 bottom right sub matrices of A and B. Thus:
0452         //  = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1)
0453         //  = AB^-1AB^-1 + det(M) - tr(M)(AB^-1)
0454         // Since we are only interested in having x, y, z with a correct ratio, we have:
0455         const Scalar
0456           a11 = m_S.coeff(f,f),     a12 = m_S.coeff(f,f+1),
0457           a21 = m_S.coeff(f+1,f),   a22 = m_S.coeff(f+1,f+1),
0458                                     a32 = m_S.coeff(f+2,f+1),
0459 
0460           a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l),
0461           a98 = m_S.coeff(l,l-1),   a99 = m_S.coeff(l,l),
0462 
0463           b11 = m_T.coeff(f,f),     b12 = m_T.coeff(f,f+1),
0464                                     b22 = m_T.coeff(f+1,f+1),
0465 
0466           b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l),
0467                                     b99 = m_T.coeff(l,l);
0468 
0469         x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21)
0470           + a12/b22 - (a11/b11)*(b12/b22);
0471         y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99);
0472         z = a32/b22;
0473       }
0474 
0475       JRs G;
0476 
0477       for (Index k=f; k<=l-2; k++)
0478       {
0479         // variables for Householder reflections
0480         Vector2s essential2;
0481         Scalar tau, beta;
0482 
0483         Vector3s hr(x,y,z);
0484 
0485         // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1)
0486         hr.makeHouseholderInPlace(tau, beta);
0487         essential2 = hr.template bottomRows<2>();
0488         Index fc=(std::max)(k-1,Index(0));  // first col to update
0489         m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
0490         m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
0491         if (m_computeQZ)
0492           m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data());
0493         if (k>f)
0494           m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0);
0495 
0496         // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k)
0497         hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1);
0498         hr.makeHouseholderInPlace(tau, beta);
0499         essential2 = hr.template bottomRows<2>();
0500         {
0501           Index lr = (std::min)(k+4,dim); // last row to update
0502           Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr);
0503           // S
0504           tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2;
0505           tmp += m_S.col(k+2).head(lr);
0506           m_S.col(k+2).head(lr) -= tau*tmp;
0507           m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
0508           // T
0509           tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2;
0510           tmp += m_T.col(k+2).head(lr);
0511           m_T.col(k+2).head(lr) -= tau*tmp;
0512           m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
0513         }
0514         if (m_computeQZ)
0515         {
0516           // Z
0517           Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim);
0518           tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k));
0519           tmp += m_Z.row(k+2);
0520           m_Z.row(k+2) -= tau*tmp;
0521           m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp);
0522         }
0523         m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0);
0524 
0525         // Z_{k2} to annihilate T(k+1,k)
0526         G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k));
0527         m_S.applyOnTheRight(k+1,k,G);
0528         m_T.applyOnTheRight(k+1,k,G);
0529         // update Z
0530         if (m_computeQZ)
0531           m_Z.applyOnTheLeft(k+1,k,G.adjoint());
0532         m_T.coeffRef(k+1,k) = Scalar(0.0);
0533 
0534         // update x,y,z
0535         x = m_S.coeff(k+1,k);
0536         y = m_S.coeff(k+2,k);
0537         if (k < l-2)
0538           z = m_S.coeff(k+3,k);
0539       } // loop over k
0540 
0541       // Q_{n-1} to annihilate y = S(l,l-2)
0542       G.makeGivens(x,y);
0543       m_S.applyOnTheLeft(l-1,l,G.adjoint());
0544       m_T.applyOnTheLeft(l-1,l,G.adjoint());
0545       if (m_computeQZ)
0546         m_Q.applyOnTheRight(l-1,l,G);
0547       m_S.coeffRef(l,l-2) = Scalar(0.0);
0548 
0549       // Z_{n-1} to annihilate T(l,l-1)
0550       G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1));
0551       m_S.applyOnTheRight(l,l-1,G);
0552       m_T.applyOnTheRight(l,l-1,G);
0553       if (m_computeQZ)
0554         m_Z.applyOnTheLeft(l,l-1,G.adjoint());
0555       m_T.coeffRef(l,l-1) = Scalar(0.0);
0556     }
0557 
0558   template<typename MatrixType>
0559     RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ)
0560     {
0561 
0562       const Index dim = A_in.cols();
0563 
0564       eigen_assert (A_in.rows()==dim && A_in.cols()==dim 
0565           && B_in.rows()==dim && B_in.cols()==dim 
0566           && "Need square matrices of the same dimension");
0567 
0568       m_isInitialized = true;
0569       m_computeQZ = computeQZ;
0570       m_S = A_in; m_T = B_in;
0571       m_workspace.resize(dim*2);
0572       m_global_iter = 0;
0573 
0574       // entrance point: hessenberg triangular decomposition
0575       hessenbergTriangular();
0576       // compute L1 vector norms of T, S into m_normOfS, m_normOfT
0577       computeNorms();
0578 
0579       Index l = dim-1, 
0580             f, 
0581             local_iter = 0;
0582 
0583       while (l>0 && local_iter<m_maxIters)
0584       {
0585         f = findSmallSubdiagEntry(l);
0586         // now rows and columns f..l (including) decouple from the rest of the problem
0587         if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0);
0588         if (f == l) // One root found
0589         {
0590           l--;
0591           local_iter = 0;
0592         }
0593         else if (f == l-1) // Two roots found
0594         {
0595           splitOffTwoRows(f);
0596           l -= 2;
0597           local_iter = 0;
0598         }
0599         else // No convergence yet
0600         {
0601           // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations
0602           Index z = findSmallDiagEntry(f,l);
0603           if (z>=f)
0604           {
0605             // zero found
0606             pushDownZero(z,f,l);
0607           }
0608           else
0609           {
0610             // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg 
0611             // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to
0612             // apply a QR-like iteration to rows and columns f..l.
0613             step(f,l, local_iter);
0614             local_iter++;
0615             m_global_iter++;
0616           }
0617         }
0618       }
0619       // check if we converged before reaching iterations limit
0620       m_info = (local_iter<m_maxIters) ? Success : NoConvergence;
0621 
0622       // For each non triangular 2x2 diagonal block of S,
0623       //    reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD.
0624       // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors,
0625       // and is in par with Lapack/Matlab QZ.
0626       if(m_info==Success)
0627       {
0628         for(Index i=0; i<dim-1; ++i)
0629         {
0630           if(m_S.coeff(i+1, i) != Scalar(0))
0631           {
0632             JacobiRotation<Scalar> j_left, j_right;
0633             internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);
0634 
0635             // Apply resulting Jacobi rotations
0636             m_S.applyOnTheLeft(i,i+1,j_left);
0637             m_S.applyOnTheRight(i,i+1,j_right);
0638             m_T.applyOnTheLeft(i,i+1,j_left);
0639             m_T.applyOnTheRight(i,i+1,j_right);
0640             m_T(i+1,i) = m_T(i,i+1) = Scalar(0);
0641 
0642             if(m_computeQZ) {
0643               m_Q.applyOnTheRight(i,i+1,j_left.transpose());
0644               m_Z.applyOnTheLeft(i,i+1,j_right.transpose());
0645             }
0646 
0647             i++;
0648           }
0649         }
0650       }
0651 
0652       return *this;
0653     } // end compute
0654 
0655 } // end namespace Eigen
0656 
0657 #endif //EIGEN_REAL_QZ
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!