diff options
Diffstat (limited to 'Eigen/src/Eigenvalues/RealQZ.h')
| -rw-r--r-- | Eigen/src/Eigenvalues/RealQZ.h | 48 |
1 files changed, 39 insertions, 9 deletions
diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h index fba6f1d77..b3a910dd9 100644 --- a/Eigen/src/Eigenvalues/RealQZ.h +++ b/Eigen/src/Eigenvalues/RealQZ.h @@ -67,7 +67,7 @@ namespace Eigen { }; typedef typename MatrixType::Scalar Scalar; typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar; - typedef typename MatrixType::Index Index; + 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; @@ -83,7 +83,7 @@ namespace Eigen { * * \sa compute() for an example. */ - RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : + explicit RealQZ(Index size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime) : m_S(size, size), m_T(size, size), m_Q(size, size), @@ -244,7 +244,7 @@ namespace Eigen { if (m_computeQZ) m_Q.applyOnTheRight(i-1,i,G); } - // update Q + // 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)); @@ -276,7 +276,7 @@ namespace Eigen { /** \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) + inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) { using std::abs; Index res = iu; @@ -294,7 +294,7 @@ namespace Eigen { /** \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) + inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) { using std::abs; Index res = l; @@ -315,8 +315,8 @@ namespace Eigen { const Index dim=m_S.cols(); if (abs(m_S.coeff(i+1,i))==Scalar(0)) return; - Index z = findSmallDiagEntry(i,i+1); - if (z==i-1) + 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>(). @@ -352,7 +352,7 @@ namespace Eigen { } else { - pushDownZero(z,i,i+1); + pushDownZero(j,i,i+1); } } @@ -552,7 +552,6 @@ namespace Eigen { 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) { @@ -616,6 +615,37 @@ namespace Eigen { } // 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 |