11 #ifndef EIGEN_REAL_SCHUR_H
12 #define EIGEN_REAL_SCHUR_H
14 #include "./HessenbergDecomposition.h"
57 typedef _MatrixType MatrixType;
59 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
60 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
61 Options = MatrixType::Options,
62 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
65 typedef typename MatrixType::Scalar Scalar;
66 typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
67 typedef typename MatrixType::Index Index;
86 m_workspaceVector(size),
88 m_isInitialized(false),
89 m_matUisUptodate(false)
102 RealSchur(
const MatrixType& matrix,
bool computeU =
true)
103 : m_matT(matrix.rows(),matrix.cols()),
104 m_matU(matrix.rows(),matrix.cols()),
105 m_workspaceVector(matrix.rows()),
106 m_hess(matrix.rows()),
107 m_isInitialized(false),
108 m_matUisUptodate(false)
126 eigen_assert(m_isInitialized &&
"RealSchur is not initialized.");
127 eigen_assert(m_matUisUptodate &&
"The matrix U has not been computed during the RealSchur decomposition.");
143 eigen_assert(m_isInitialized &&
"RealSchur is not initialized.");
172 eigen_assert(m_isInitialized &&
"RealSchur is not initialized.");
189 bool m_isInitialized;
190 bool m_matUisUptodate;
194 Scalar computeNormOfT();
195 Index findSmallSubdiagEntry(Index iu, Scalar norm);
196 void splitOffTwoRows(Index iu,
bool computeU, Scalar exshift);
197 void computeShift(Index iu, Index iter, Scalar& exshift,
Vector3s& shiftInfo);
198 void initFrancisQRStep(Index il, Index iu,
const Vector3s& shiftInfo, Index& im,
Vector3s& firstHouseholderVector);
199 void performFrancisQRStep(Index il, Index im, Index iu,
bool computeU,
const Vector3s& firstHouseholderVector, Scalar* workspace);
203 template<
typename MatrixType>
206 assert(matrix.cols() == matrix.rows());
210 m_matT = m_hess.matrixH();
212 m_matU = m_hess.matrixQ();
215 m_workspaceVector.resize(m_matT.cols());
216 Scalar* workspace = &m_workspaceVector.coeffRef(0);
222 Index iu = m_matT.cols() - 1;
226 Scalar norm = computeNormOfT();
232 Index il = findSmallSubdiagEntry(iu, norm);
237 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
239 m_matT.coeffRef(iu, iu-1) = Scalar(0);
245 splitOffTwoRows(iu, computeU, exshift);
252 Vector3s firstHouseholderVector(0,0,0), shiftInfo;
253 computeShift(iu, iter, exshift, shiftInfo);
255 totalIter = totalIter + 1;
256 if (totalIter > m_maxIterations * matrix.cols())
break;
258 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
259 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
263 if(totalIter <= m_maxIterations * matrix.cols())
268 m_isInitialized =
true;
269 m_matUisUptodate = computeU;
274 template<
typename MatrixType>
277 const Index size = m_matT.cols();
282 for (Index j = 0; j < size; ++j)
283 norm += m_matT.row(j).segment((std::max)(j-1,Index(0)), size-(std::max)(j-1,Index(0))).cwiseAbs().sum();
288 template<
typename MatrixType>
289 inline typename MatrixType::Index RealSchur<MatrixType>::findSmallSubdiagEntry(Index iu, Scalar norm)
294 Scalar s = internal::abs(m_matT.coeff(res-1,res-1)) + internal::abs(m_matT.coeff(res,res));
297 if (internal::abs(m_matT.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
305 template<
typename MatrixType>
306 inline void RealSchur<MatrixType>::splitOffTwoRows(Index iu,
bool computeU, Scalar exshift)
308 const Index size = m_matT.cols();
312 Scalar p = Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
313 Scalar q = p * p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
314 m_matT.coeffRef(iu,iu) += exshift;
315 m_matT.coeffRef(iu-1,iu-1) += exshift;
319 Scalar z = internal::sqrt(internal::abs(q));
320 JacobiRotation<Scalar> rot;
322 rot.makeGivens(p + z, m_matT.coeff(iu, iu-1));
324 rot.makeGivens(p - z, m_matT.coeff(iu, iu-1));
326 m_matT.rightCols(size-iu+1).applyOnTheLeft(iu-1, iu, rot.adjoint());
327 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu, rot);
328 m_matT.coeffRef(iu, iu-1) = Scalar(0);
330 m_matU.applyOnTheRight(iu-1, iu, rot);
334 m_matT.coeffRef(iu-1, iu-2) = Scalar(0);
338 template<
typename MatrixType>
339 inline void RealSchur<MatrixType>::computeShift(Index iu, Index iter, Scalar& exshift, Vector3s& shiftInfo)
341 shiftInfo.coeffRef(0) = m_matT.coeff(iu,iu);
342 shiftInfo.coeffRef(1) = m_matT.coeff(iu-1,iu-1);
343 shiftInfo.coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
348 exshift += shiftInfo.coeff(0);
349 for (Index i = 0; i <= iu; ++i)
350 m_matT.coeffRef(i,i) -= shiftInfo.coeff(0);
351 Scalar s = internal::abs(m_matT.coeff(iu,iu-1)) + internal::abs(m_matT.coeff(iu-1,iu-2));
352 shiftInfo.coeffRef(0) = Scalar(0.75) * s;
353 shiftInfo.coeffRef(1) = Scalar(0.75) * s;
354 shiftInfo.coeffRef(2) = Scalar(-0.4375) * s * s;
360 Scalar s = (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
361 s = s * s + shiftInfo.coeff(2);
364 s = internal::sqrt(s);
365 if (shiftInfo.coeff(1) < shiftInfo.coeff(0))
367 s = s + (shiftInfo.coeff(1) - shiftInfo.coeff(0)) / Scalar(2.0);
368 s = shiftInfo.coeff(0) - shiftInfo.coeff(2) / s;
370 for (Index i = 0; i <= iu; ++i)
371 m_matT.coeffRef(i,i) -= s;
372 shiftInfo.setConstant(Scalar(0.964));
378 template<
typename MatrixType>
379 inline void RealSchur<MatrixType>::initFrancisQRStep(Index il, Index iu,
const Vector3s& shiftInfo, Index& im, Vector3s& firstHouseholderVector)
381 Vector3s& v = firstHouseholderVector;
383 for (im = iu-2; im >= il; --im)
385 const Scalar Tmm = m_matT.coeff(im,im);
386 const Scalar r = shiftInfo.coeff(0) - Tmm;
387 const Scalar s = shiftInfo.coeff(1) - Tmm;
388 v.coeffRef(0) = (r * s - shiftInfo.coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
389 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r - s;
390 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
394 const Scalar lhs = m_matT.coeff(im,im-1) * (internal::abs(v.coeff(1)) + internal::abs(v.coeff(2)));
395 const Scalar rhs = v.coeff(0) * (internal::abs(m_matT.coeff(im-1,im-1)) + internal::abs(Tmm) + internal::abs(m_matT.coeff(im+1,im+1)));
396 if (internal::abs(lhs) < NumTraits<Scalar>::epsilon() * rhs)
404 template<
typename MatrixType>
405 inline void RealSchur<MatrixType>::performFrancisQRStep(Index il, Index im, Index iu,
bool computeU,
const Vector3s& firstHouseholderVector, Scalar* workspace)
410 const Index size = m_matT.cols();
412 for (Index k = im; k <= iu-2; ++k)
414 bool firstIteration = (k == im);
418 v = firstHouseholderVector;
420 v = m_matT.template block<3,1>(k,k-1);
423 Matrix<Scalar, 2, 1> ess;
424 v.makeHouseholder(ess, tau, beta);
426 if (beta != Scalar(0))
428 if (firstIteration && k > il)
429 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
430 else if (!firstIteration)
431 m_matT.coeffRef(k,k-1) = beta;
434 m_matT.block(k, k, 3, size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
435 m_matT.block(0, k, (std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
437 m_matU.block(0, k, size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
441 Matrix<Scalar, 2, 1> v = m_matT.template block<2,1>(iu-1, iu-2);
443 Matrix<Scalar, 1, 1> ess;
444 v.makeHouseholder(ess, tau, beta);
446 if (beta != Scalar(0))
448 m_matT.coeffRef(iu-1, iu-2) = beta;
449 m_matT.block(iu-1, iu-1, 2, size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
450 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
452 m_matU.block(0, iu-1, size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
456 for (Index i = im+2; i <= iu; ++i)
458 m_matT.coeffRef(i,i-2) = Scalar(0);
460 m_matT.coeffRef(i,i-3) = Scalar(0);
466 #endif // EIGEN_REAL_SCHUR_H