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