Commits

Alexey Korepanov  committed a3a2629

+Gael comments

  • Participants
  • Parent commits 7a28d80

Comments (0)

Files changed (1)

File Eigen/src/Eigenvalues/RealQZ.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.
-  *
-  * I should add an example of usage of this class, but I don't exactly know
-  * how.
-  *
-  * \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
-  */
+  /** \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.
+   *
+   * I should add an example of usage of this class, but I don't exactly know
+   * how.
+   *
+   * \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 typename MatrixType::Index Index;
+  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 typename MatrixType::Index Index;
 
-            typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
-            typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;
+      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.
-             */
-            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_isInitialized(false)
-                { }
+      /** \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.
+       */
+      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_isInitialized(false)
+        { }
 
-            /** \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_isInitialized(false) {
-                    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 succesful, \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;
-            }
-
-            /** \brief Maximum number of iterations.
-             *
-             * Maximum number of iterations allowed for an eigenvalue to converge. 
-             */
-            static const Index m_max_iter = 400;
-
-        private:
-
-            MatrixType m_S, m_T, m_Q, m_Z;
-            ColumnVectorType m_workspace;
-            ComputationInfo m_info;
-            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)
-                    G.makeGivens(m_S(i-1, j), m_S(i,j));
-                    m_S.applyOnTheLeft(i-1,i,G.adjoint());
-                    m_T.applyOnTheLeft(i-1,i,G.adjoint());
-                    m_S(i,j) = 0.0;
-                    // update Q
-                    if (m_computeQZ)
-                        m_Q.applyOnTheRight(i-1,i,G);
-                    // kill T(i,i-1)
-                    G.makeGivens(m_T(i,i), m_T(i,i-1));
-                    m_S.applyOnTheRight(i,i-1,G);
-                    m_T.applyOnTheRight(i,i-1,G);
-                    m_T(i,i-1) = 0.0;
-                    // update Z
-                    if (m_computeQZ)
-                        m_Z.applyOnTheLeft(i,i-1,G.adjoint());
-                }
-            }
+      /** \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_isInitialized(false) {
+          compute(A, B, computeQZ);
         }
 
-    /** \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);
-            m_normOfT = Scalar(0);
-            for (Index j = 0; j < size; ++j) {
-                Index row_start = (std::max)(j-1,Index(0));
-                m_normOfS += m_S.row(j).segment(row_start, size - row_start).cwiseAbs().sum();
-                m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
-            }
+      /** \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 succesful, \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;
+      }
+
+      /** \brief Maximum number of iterations.
+       *
+       * Maximum number of iterations allowed for an eigenvalue to converge. 
+       */
+      static const Index m_max_iter = 400;
+
+    private:
+
+      MatrixType m_S, m_T, m_Q, m_Z;
+      ColumnVectorType m_workspace;
+      ComputationInfo m_info;
+      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)
+          G.makeGivens(m_S.coeff(i-1, j), m_S.coeff(i,j));
+          m_S.applyOnTheLeft(i-1,i,G.adjoint());
+          m_T.applyOnTheLeft(i-1,i,G.adjoint());
+          m_S.coeffRef(i,j) = Scalar(0.0);
+          // update Q
+          if (m_computeQZ)
+            m_Q.applyOnTheRight(i-1,i,G);
+          // kill T(i,i-1)
+          G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1));
+          m_S.applyOnTheRight(i,i-1,G);
+          m_T.applyOnTheRight(i,i-1,G);
+          m_T.coeffRef(i,i-1) = Scalar(0.0);
+          // 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) {
+        Index row_start = (std::max)(j-1,Index(0));
+        m_normOfS += m_S.row(j).segment(row_start, size - row_start).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 typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) {
+      Index res = iu;
+      while (res > 0) {
+        Scalar s = internal::abs(m_S.coeff(res-1,res-1)) + internal::abs(m_S.coeff(res,res));
+        if (s == Scalar(0.0))
+          s = m_normOfS;
+        if (internal::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 typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) {
+      Index res = l;
+      while (res >= f) {
+        if (internal::abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
+          break;
+        res--;
+      }
+      return res;
+    }
+
+  /** \internal decouple 2x2 diagonal block in rows iu, iu+1 if eigenvalues are real */
+  template<typename MatrixType>
+    inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) {
+      const Index dim=m_S.cols();
+      if (internal::abs(m_S.coeff(i+1,1)==Scalar(0)))
+        return;
+      Index z = findSmallDiagEntry(i,i+1);
+      if (z==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 = internal::sqrt(q);
+          // QR for ABi - lambda I
+          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(z,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.bottomRows(dim-zz).applyOnTheRight(zz, zz-1,G);
+          m_T.bottomRows(dim-zz).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 */
+  template<typename MatrixType>
+    inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) {
+      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 = internal::abs(a87*b77i) + internal::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 {
+        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), b11i=1.0/b11, b12=m_T.coeff(f,f+1),
+          b22i=Scalar(1.0)/m_T.coeff(f+1,f+1),
+          b88i=Scalar(1.0)/m_T.coeff(l-1,l-1), b89=m_T.coeff(l-1,l),
+          b99i=Scalar(1.0)/m_T.coeff(l,l);
+        x = ( (a88*b88i - a11*b11i)*(a99*b99i - a11*b11i) - (a89*b99i)*(a98*b88i) + (a98*b88i)*(b89*b99i)*(a11*b11i) ) * (b11/a21)
+          + a12*b22i - (a11*b11i)*(b12*b22i);
+        y = (a22*b22i-a11*b11i) - (a21*b11i)*(b12*b22i) - (a88*b88i-a11*b11i) - (a99*b99i-a11*b11i) + (a98*b88i)*(b89*b99i);
+        z = a32*b22i;
+      }
+
+      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
+        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+1,k-1) = Scalar(0.0);
+          m_S.coeffRef(k+2,k-1) = Scalar(0.0);
         }
 
-    
-    /** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
-    template<typename MatrixType>
-        inline typename MatrixType::Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) {
-            Index res = iu;
-            while (res > 0) {
-                Scalar s = internal::abs(m_S(res-1,res-1)) + internal::abs(m_S(res,res));
-                if (s == 0.0)
-                    s = m_normOfS;
-                if (internal::abs(m_S(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
-                    break;
-                res--;
-            }
-            return res;
+        // Z_{k1}
+        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();
         }
-    
-    /** \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 typename MatrixType::Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) {
-            Index res = l;
-            while (res >= f) {
-                if (internal::abs(m_T(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
-                    break;
-                res--;
-            }
-            return res;
+        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) = Scalar(0.0);
+        m_T.coeffRef(k+2,k+1) = Scalar(0.0);
 
-    /** \internal decouple 2x2 diagonal block in rows iu, iu+1 if eigenvalues are real */
-    template<typename MatrixType>
-        inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) {
-            if (internal::abs(m_S(i+1,1)==Scalar(0)))
-                return;
-            Index z = findSmallDiagEntry(i,i+1);
-            if (z==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 = 0.5*(STi(0,0)-STi(1,1));
-                Scalar q = p*p + STi(1,0)*STi(0,1);
-                if (q>=0) {
-                    Scalar z = internal::sqrt(q);
-                    // QR for ABi - lambda I
-                    JRs G;
-                    if (p>=0)
-                        G.makeGivens(p + z, STi(1,0));
-                    else
-                        G.makeGivens(p - z, STi(1,0));
-                    m_S.applyOnTheLeft(i,i+1,G.adjoint());
-                    m_T.applyOnTheLeft(i,i+1,G.adjoint());
-                    // update Q
-                    if (m_computeQZ)
-                        m_Q.applyOnTheRight(i,i+1,G);
+        // Z_{k2}
+        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);
 
-                    G.makeGivens(m_T(i+1,i+1), m_T(i+1,i));
-                    m_S.applyOnTheRight(i+1,i,G);
-                    m_T.applyOnTheRight(i+1,i,G);
-                    // update Z
-                    if (m_computeQZ)
-                        m_Z.applyOnTheLeft(i+1,i,G.adjoint());
+        // 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
 
-                    m_S(i+1,i) = 0.0;
-                    m_T(i+1,i) = 0.0;
-                }
-            } else {
-                pushDownZero(z,i,i+1);
-            }
+      // Q_{n-1}
+      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}
+      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();
+
+      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_max_iter) {
+        f = findSmallSubdiagEntry(l);
+        if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0);
+        if (f == l) {
+          l --;
+          local_iter = 0;
+        } else if (f == l-1) {
+          splitOffTwoRows(f);
+          l -= 2;
+          local_iter = 0;
+        } else {
+          Index z = findSmallDiagEntry(f,l);
+          if (z>=f) {
+            // zero found
+            pushDownZero(z,f,l);
+          } else {
+            // QR-like iteration
+            step(f,l, local_iter);
+            local_iter++;
+            m_global_iter++;
+          }
         }
-    
-    /** \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(zz, zz+1), m_T(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(zz+1,zz+1) = 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(zz+1, zz), m_S(zz+1,zz-1));
-                    m_S.bottomRows(dim-zz).applyOnTheRight(zz, zz-1,G);
-                    m_T.bottomRows(dim-zz).applyOnTheRight(zz, zz-1,G);
-                    m_S(zz+1,zz-1) = 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(l,l), m_S(l,l-1));
-            m_S.applyOnTheRight(l,l-1,G);
-            m_T.applyOnTheRight(l,l-1,G);
-            m_S(l,l-1)=0.0;
-            // update Z
-            if (m_computeQZ)
-                m_Z.applyOnTheLeft(l,l-1,G.adjoint());
-        }
-    
-    /** \internal QR-like iterative step */
-    template<typename MatrixType>
-        inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) {
-            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(f+0,f+0), a12=m_S(f+0,f+1),
-                    a21=m_S(f+1,f+0), a22=m_S(f+1,f+1), a32=m_S(f+2,f+1),
-                    b12=m_T(f+0,f+1),
-                    b11i=1.0/m_T(f+0,f+0),
-                    b22i=1.0/m_T(f+1,f+1),
-                    a87=m_S(l-1,l-2),
-                    a98=m_S(l-0,l-1),
-                    b77i=1.0/m_T(l-2,l-2),
-                    b88i=1.0/m_T(l-1,l-1);
-                Scalar ss = internal::abs(a87*b77i) + internal::abs(a98*b88i),
-                       lpl = 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==20 || iter==23 || iter == 27) {
-                // another exceptional shift
-                x = m_S(f,f)/m_T(f,f)-m_S(l,l)/m_T(l,l) + m_S(l,l-1)*m_T(l-1,l)/(m_T(l-1,l-1)*m_T(l,l));
-                y = m_S(f+1,f)/m_T(f,f);
-                z = 0;
-            } else if (iter>31 && !(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 {
-                const Scalar 
-                    a11=m_S(f,f), a12=m_S(f,f+1), 
-                    a21=m_S(f+1,f), a22=m_S(f+1,f+1),
-                    a32=m_S(f+2,f+1),
-                    a88=m_S(l-1,l-1), a89=m_S(l-1,l),
-                    a98=m_S(l,l-1), a99=m_S(l,l),
-                    b11=m_T(f,f), b11i=1.0/b11, b12=m_T(f,f+1),
-                    b22i=1.0/m_T(f+1,f+1),
-                    b88i=1.0/m_T(l-1,l-1), b89=m_T(l-1,l),
-                    b99i=1.0/m_T(l,l);
-                x = ( (a88*b88i - a11*b11i)*(a99*b99i - a11*b11i) - (a89*b99i)*(a98*b88i) + (a98*b88i)*(b89*b99i)*(a11*b11i) ) * (b11/a21)
-                    + a12*b22i - (a11*b11i)*(b12*b22i);
-                y = (a22*b22i-a11*b11i) - (a21*b11i)*(b12*b22i) - (a88*b88i-a11*b11i) - (a99*b99i-a11*b11i) + (a98*b88i)*(b89*b99i);
-                z = a32*b22i;
-            }
-
-            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
-                hr.makeHouseholderInPlace(tau, beta);
-                essential2 = hr.template bottomRows<2>();
-                Index fc=std::max(k-1,Index(0));  // first col to update
-                m_S.block(k,fc,3,dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
-                m_T.block(k,fc,3,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(k+1,k-1) = 0.0;
-                    m_S(k+2,k-1) = 0.0;
-                }
-
-                // Z_{k1}
-                hr << m_T(k+2,k+2),m_T(k+2,k),m_T(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.block(0,k,lr,2) * essential2;
-                    tmp += m_S.block(0,k+2,lr,1);
-                    m_S.block(0,k+2,lr,1) -= tau*tmp;
-                    m_S.block(0,k,lr,2) -= (tau*tmp) * essential2.adjoint();
-                    // T
-                    tmp = m_T.block(0,k,lr,2) * essential2;
-                    tmp += m_T.block(0,k+2,lr,1);
-                    m_T.block(0,k+2,lr,1) -= tau*tmp;
-                    m_T.block(0,k,lr,2) -= (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(k+2,k) = 0.0;
-                m_T(k+2,k+1) = 0.0;
-
-                // Z_{k2}
-                G.makeGivens(m_T(k+1,k+1), m_T(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(k+1,k) = 0.0;
-
-                // update x,y,z
-                x = m_S(k+1,k);
-                y = m_S(k+2,k);
-                if (k < l-2)
-                    z = m_S(k+3,k);
-            } // loop over k
-
-            // Q_{n-1}
-            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(l,l-2) = 0.0;
-
-            // Z_{n-1}
-            G.makeGivens(m_T(l,l),m_T(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(l,l-1) = 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();
-            
-            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_max_iter) {
-                f = findSmallSubdiagEntry(l);
-                if (f>0) m_S(f,f-1) = 0.0;
-                if (f == l) {
-                    l --;
-                    local_iter = 0;
-                } else if (f == l-1) {
-                    splitOffTwoRows(f);
-                    l -= 2;
-                    local_iter = 0;
-                } else {
-                    Index z = findSmallDiagEntry(f,l);
-                    if (z>=f) {
-                        // zero found
-                        pushDownZero(z,f,l);
-                    } else {
-                        // QR-like iteration
-                        step(f,l, local_iter);
-                        local_iter++;
-                        m_global_iter++;
-                    }
-                }
-            }
-            // check if we converged before reaching iterations limit
-            if (local_iter<m_max_iter) {
-                m_info = Success;
-            } else {
-                m_info = NoConvergence;
-            }
-            return *this;
-        } // end compute
+      }
+      // check if we converged before reaching iterations limit
+      if (local_iter<m_max_iter) {
+        m_info = Success;
+      } else {
+        m_info = NoConvergence;
+      }
+      return *this;
+    } // end compute
 
 
 } // end namespace Eigen