aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/Eigenvalues
diff options
context:
space:
mode:
authorGravatar Gael Guennebaud <g.gael@free.fr>2016-06-09 17:11:03 +0200
committerGravatar Gael Guennebaud <g.gael@free.fr>2016-06-09 17:11:03 +0200
commitc1f9ca925405c8fad126f327b4bdca7f983fc96e (patch)
tree06f6cb2be66d89cc553426d335a8c3ee6b0c4689 /Eigen/src/Eigenvalues
parent15890c304edbccedc8a989468ed3fc475f428059 (diff)
Update RealQZ to reduce 2x2 diagonal block of T corresponding to non reduced diagonal block of S to positive diagonal form.
This step involve a real 2x2 SVD problem. The respective routine is thus in src/misc/ to be shared by both EVD and AVD modules.
Diffstat (limited to 'Eigen/src/Eigenvalues')
-rw-r--r--Eigen/src/Eigenvalues/RealQZ.h32
1 files changed, 31 insertions, 1 deletions
diff --git a/Eigen/src/Eigenvalues/RealQZ.h b/Eigen/src/Eigenvalues/RealQZ.h
index a62071d42..c4715b954 100644
--- a/Eigen/src/Eigenvalues/RealQZ.h
+++ b/Eigen/src/Eigenvalues/RealQZ.h
@@ -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_T.applyOnTheLeft(i,i+1,j_left);
+ m_T.applyOnTheRight(i,i+1,j_right);
+ m_S.applyOnTheLeft(i,i+1,j_left);
+ m_S.applyOnTheRight(i,i+1,j_right);
+ 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