10 #ifndef EIGEN_REAL_QZ_H
11 #define EIGEN_REAL_QZ_H
57 template<
typename _MatrixType>
class RealQZ
219 template<
typename MatrixType>
223 const Index dim = m_S.cols();
228 m_T.template triangularView<StrictlyLower>().setZero();
231 m_S.applyOnTheLeft(m_Q.adjoint());
234 m_Z = MatrixType::Identity(dim,dim);
240 if(m_S.coeff(
i,
j) != 0)
242 G.makeGivens(m_S.coeff(
i-1,
j), m_S.coeff(
i,
j), &m_S.coeffRef(
i-1,
j));
244 m_S.rightCols(dim-
j-1).applyOnTheLeft(
i-1,
i,
G.adjoint());
245 m_T.rightCols(dim-
i+1).applyOnTheLeft(
i-1,
i,
G.adjoint());
248 m_Q.applyOnTheRight(
i-1,
i,
G);
253 G.makeGivens(m_T.coeff(
i,
i), m_T.coeff(
i,
i-1), &m_T.coeffRef(
i,
i));
255 m_S.applyOnTheRight(
i,
i-1,
G);
256 m_T.topRows(
i).applyOnTheRight(
i,
i-1,
G);
259 m_Z.applyOnTheLeft(
i,
i-1,
G.adjoint());
266 template<
typename MatrixType>
274 m_normOfS += m_S.col(
j).segment(0, (
std::min)(
size,
j+2)).cwiseAbs().sum();
275 m_normOfT += m_T.row(
j).segment(
j,
size -
j).cwiseAbs().sum();
281 template<
typename MatrixType>
299 template<
typename MatrixType>
313 template<
typename MatrixType>
318 const Index dim=m_S.cols();
321 Index j = findSmallDiagEntry(
i,
i+1);
325 Matrix2s STi = m_T.template block<2,2>(
i,
i).template triangularView<Upper>().
326 template solve<OnTheRight>(m_S.template block<2,2>(
i,
i));
336 G.makeGivens(
p +
z, STi(1,0));
338 G.makeGivens(
p -
z, STi(1,0));
339 m_S.rightCols(dim-
i).applyOnTheLeft(
i,
i+1,
G.adjoint());
340 m_T.rightCols(dim-
i).applyOnTheLeft(
i,
i+1,
G.adjoint());
343 m_Q.applyOnTheRight(
i,
i+1,
G);
345 G.makeGivens(m_T.coeff(
i+1,
i+1), m_T.coeff(
i+1,
i));
346 m_S.topRows(
i+2).applyOnTheRight(
i+1,
i,
G);
347 m_T.topRows(
i+2).applyOnTheRight(
i+1,
i,
G);
350 m_Z.applyOnTheLeft(
i+1,
i,
G.adjoint());
358 pushDownZero(
j,
i,
i+1);
363 template<
typename MatrixType>
367 const Index dim = m_S.cols();
371 Index firstColS = zz>
f ? (zz-1) : zz;
372 G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1));
373 m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,
G.adjoint());
374 m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,
G.adjoint());
375 m_T.coeffRef(zz+1,zz+1) =
Scalar(0.0);
378 m_Q.applyOnTheRight(zz,zz+1,
G);
382 G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1));
383 m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,
G);
384 m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,
G);
385 m_S.coeffRef(zz+1,zz-1) =
Scalar(0.0);
388 m_Z.applyOnTheLeft(zz,zz-1,
G.adjoint());
392 G.makeGivens(m_S.coeff(
l,
l), m_S.coeff(
l,
l-1));
393 m_S.applyOnTheRight(
l,
l-1,
G);
394 m_T.applyOnTheRight(
l,
l-1,
G);
398 m_Z.applyOnTheLeft(
l,
l-1,
G.adjoint());
402 template<
typename MatrixType>
406 const Index dim = m_S.cols();
414 a11=m_S.coeff(
f+0,
f+0), a12=m_S.coeff(
f+0,
f+1),
415 a21=m_S.coeff(
f+1,
f+0), a22=m_S.coeff(
f+1,
f+1), a32=m_S.coeff(
f+2,
f+1),
416 b12=m_T.coeff(
f+0,
f+1),
417 b11i=
Scalar(1.0)/m_T.coeff(
f+0,
f+0),
418 b22i=
Scalar(1.0)/m_T.coeff(
f+1,
f+1),
419 a87=m_S.coeff(
l-1,
l-2),
420 a98=m_S.coeff(
l-0,
l-1),
421 b77i=
Scalar(1.0)/m_T.coeff(
l-2,
l-2),
422 b88i=
Scalar(1.0)/m_T.coeff(
l-1,
l-1);
426 x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i
427 - a11*a21*b12*b11i*b11i*b22i;
428 y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i
429 - a21*a21*b12*b11i*b11i*b22i;
430 z = a21*a32*b11i*b22i;
435 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) /
436 (m_T.coeff(
l-1,
l-1)*m_T.coeff(
l,
l));
437 y = m_S.coeff(
f+1,
f)/m_T.coeff(
f,
f);
443 x = internal::random<Scalar>(-1.0,1.0);
444 y = internal::random<Scalar>(-1.0,1.0);
445 z = internal::random<Scalar>(-1.0,1.0);
456 a11 = m_S.coeff(
f,
f), a12 = m_S.coeff(
f,
f+1),
457 a21 = m_S.coeff(
f+1,
f), a22 = m_S.coeff(
f+1,
f+1),
458 a32 = m_S.coeff(
f+2,
f+1),
460 a88 = m_S.coeff(
l-1,
l-1), a89 = m_S.coeff(
l-1,
l),
461 a98 = m_S.coeff(
l,
l-1), a99 = m_S.coeff(
l,
l),
463 b11 = m_T.coeff(
f,
f), b12 = m_T.coeff(
f,
f+1),
464 b22 = m_T.coeff(
f+1,
f+1),
466 b88 = m_T.coeff(
l-1,
l-1), b89 = m_T.coeff(
l-1,
l),
467 b99 = m_T.coeff(
l,
l);
469 x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21)
470 + a12/b22 - (a11/b11)*(b12/b22);
471 y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99);
486 hr.makeHouseholderInPlace(tau,
beta);
487 essential2 = hr.template bottomRows<2>();
489 m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.
data());
490 m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.
data());
492 m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.
data());
494 m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) =
Scalar(0.0);
497 hr << m_T.
coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1);
498 hr.makeHouseholderInPlace(tau,
beta);
499 essential2 = hr.template bottomRows<2>();
504 tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2;
505 tmp += m_S.col(k+2).head(lr);
506 m_S.col(k+2).head(lr) -= tau*tmp;
507 m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
509 tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2;
510 tmp += m_T.col(k+2).head(lr);
511 m_T.col(k+2).head(lr) -= tau*tmp;
512 m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
518 tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k));
520 m_Z.row(k+2) -= tau*tmp;
521 m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp);
526 G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k));
527 m_S.applyOnTheRight(k+1,k,
G);
528 m_T.applyOnTheRight(k+1,k,
G);
531 m_Z.applyOnTheLeft(k+1,k,
G.adjoint());
532 m_T.coeffRef(k+1,k) =
Scalar(0.0);
535 x = m_S.coeff(k+1,k);
536 y = m_S.coeff(k+2,k);
538 z = m_S.coeff(k+3,k);
543 m_S.applyOnTheLeft(
l-1,
l,
G.adjoint());
544 m_T.applyOnTheLeft(
l-1,
l,
G.adjoint());
546 m_Q.applyOnTheRight(
l-1,
l,
G);
550 G.makeGivens(m_T.coeff(
l,
l),m_T.coeff(
l,
l-1));
551 m_S.applyOnTheRight(
l,
l-1,
G);
552 m_T.applyOnTheRight(
l,
l-1,
G);
554 m_Z.applyOnTheLeft(
l,
l-1,
G.adjoint());
558 template<
typename MatrixType>
562 const Index dim = A_in.cols();
565 && B_in.rows()==dim && B_in.cols()==dim
566 &&
"Need square matrices of the same dimension");
568 m_isInitialized =
true;
569 m_computeQZ = computeQZ;
570 m_S = A_in; m_T = B_in;
571 m_workspace.resize(dim*2);
575 hessenbergTriangular();
583 while (
l>0 && local_iter<m_maxIters)
585 f = findSmallSubdiagEntry(
l);
587 if (
f>0) m_S.coeffRef(
f,
f-1) =
Scalar(0.0);
636 m_S.applyOnTheLeft(
i,
i+1,j_left);
637 m_S.applyOnTheRight(
i,
i+1,j_right);
638 m_T.applyOnTheLeft(
i,
i+1,j_left);
639 m_T.applyOnTheRight(
i,
i+1,j_right);
657 #endif //EIGEN_REAL_QZ