11 #ifndef EIGEN_REAL_SCHUR_H
12 #define EIGEN_REAL_SCHUR_H
103 template<
typename InputType>
169 template<
typename InputType>
189 template<
typename HessMatrixType,
typename OrthMatrixType>
247 template<
typename MatrixType>
248 template<
typename InputType>
254 Index maxIters = m_maxIters;
256 maxIters = m_maxIterationsPerRow *
matrix.rows();
259 if(
scale<considerAsZero)
265 m_isInitialized =
true;
266 m_matUisUptodate = computeU;
276 m_workspaceVector.resize(
matrix.cols());
278 m_hess.matrixQ().evalTo(m_matU, m_workspaceVector);
279 computeFromHessenberg(m_hess.matrixH(), m_matU, computeU);
285 template<
typename MatrixType>
286 template<
typename HessMatrixType,
typename OrthMatrixType>
292 m_workspaceVector.resize(m_matT.cols());
296 Index maxIters = m_maxIters;
298 maxIters = m_maxIterationsPerRow * matrixH.rows();
299 Scalar* workspace = &m_workspaceVector.coeffRef(0);
305 Index iu = m_matT.cols() - 1;
309 Scalar norm = computeNormOfT();
319 Index il = findSmallSubdiagEntry(iu,considerAsZero);
324 m_matT.coeffRef(iu,iu) = m_matT.coeff(iu,iu) + exshift;
326 m_matT.coeffRef(iu, iu-1) =
Scalar(0);
332 splitOffTwoRows(iu, computeU, exshift);
339 Vector3s firstHouseholderVector = Vector3s::Zero(), shiftInfo;
340 computeShift(iu,
iter, exshift, shiftInfo);
342 totalIter = totalIter + 1;
343 if (totalIter > maxIters)
break;
345 initFrancisQRStep(il, iu, shiftInfo, im, firstHouseholderVector);
346 performFrancisQRStep(il, im, iu, computeU, firstHouseholderVector, workspace);
350 if(totalIter <= maxIters)
355 m_isInitialized =
true;
356 m_matUisUptodate = computeU;
361 template<
typename MatrixType>
370 norm += m_matT.col(
j).segment(0, (
std::min)(
size,
j+2)).cwiseAbs().sum();
375 template<
typename MatrixType>
394 template<
typename MatrixType>
403 Scalar p =
Scalar(0.5) * (m_matT.coeff(iu-1,iu-1) - m_matT.coeff(iu,iu));
404 Scalar q =
p *
p + m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
405 m_matT.coeffRef(iu,iu) += exshift;
406 m_matT.coeffRef(iu-1,iu-1) += exshift;
413 rot.makeGivens(
p +
z, m_matT.coeff(iu, iu-1));
415 rot.makeGivens(
p -
z, m_matT.coeff(iu, iu-1));
417 m_matT.rightCols(
size-iu+1).applyOnTheLeft(iu-1, iu,
rot.adjoint());
418 m_matT.topRows(iu+1).applyOnTheRight(iu-1, iu,
rot);
419 m_matT.coeffRef(iu, iu-1) =
Scalar(0);
421 m_matU.applyOnTheRight(iu-1, iu,
rot);
425 m_matT.coeffRef(iu-1, iu-2) =
Scalar(0);
429 template<
typename MatrixType>
434 shiftInfo.
coeffRef(0) = m_matT.coeff(iu,iu);
435 shiftInfo.
coeffRef(1) = m_matT.coeff(iu-1,iu-1);
436 shiftInfo.
coeffRef(2) = m_matT.coeff(iu,iu-1) * m_matT.coeff(iu-1,iu);
441 exshift += shiftInfo.
coeff(0);
443 m_matT.coeffRef(
i,
i) -= shiftInfo.
coeff(0);
444 Scalar s =
abs(m_matT.coeff(iu,iu-1)) +
abs(m_matT.coeff(iu-1,iu-2));
464 m_matT.coeffRef(
i,
i) -=
s;
471 template<
typename MatrixType>
477 for (im = iu-2; im >= il; --im)
482 v.coeffRef(0) = (r *
s - shiftInfo.
coeff(2)) / m_matT.coeff(im+1,im) + m_matT.coeff(im,im+1);
483 v.coeffRef(1) = m_matT.coeff(im+1,im+1) - Tmm - r -
s;
484 v.coeffRef(2) = m_matT.coeff(im+2,im+1);
488 const Scalar lhs = m_matT.coeff(im,im-1) * (
abs(
v.coeff(1)) +
abs(
v.coeff(2)));
489 const Scalar rhs =
v.coeff(0) * (
abs(m_matT.coeff(im-1,im-1)) +
abs(Tmm) +
abs(m_matT.coeff(im+1,im+1)));
496 template<
typename MatrixType>
504 for (
Index k = im; k <= iu-2; ++k)
506 bool firstIteration = (k == im);
510 v = firstHouseholderVector;
512 v = m_matT.template block<3,1>(k,k-1);
516 v.makeHouseholder(ess, tau,
beta);
520 if (firstIteration && k > il)
521 m_matT.coeffRef(k,k-1) = -m_matT.coeff(k,k-1);
522 else if (!firstIteration)
523 m_matT.coeffRef(k,k-1) =
beta;
526 m_matT.block(k, k, 3,
size-k).applyHouseholderOnTheLeft(ess, tau, workspace);
527 m_matT.block(0, k, (
std::min)(iu,k+3) + 1, 3).applyHouseholderOnTheRight(ess, tau, workspace);
529 m_matU.block(0, k,
size, 3).applyHouseholderOnTheRight(ess, tau, workspace);
536 v.makeHouseholder(ess, tau,
beta);
540 m_matT.coeffRef(iu-1, iu-2) =
beta;
541 m_matT.block(iu-1, iu-1, 2,
size-iu+1).applyHouseholderOnTheLeft(ess, tau, workspace);
542 m_matT.block(0, iu-1, iu+1, 2).applyHouseholderOnTheRight(ess, tau, workspace);
544 m_matU.block(0, iu-1,
size, 2).applyHouseholderOnTheRight(ess, tau, workspace);
558 #endif // EIGEN_REAL_SCHUR_H