RealQZ.h
Go to the documentation of this file.
1 // This file is part of Eigen, a lightweight C++ template library
2 // for linear algebra.
3 //
4 // Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru>
5 //
6 // This Source Code Form is subject to the terms of the Mozilla
7 // Public License v. 2.0. If a copy of the MPL was not distributed
8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
9 
10 #ifndef EIGEN_REAL_QZ_H
11 #define EIGEN_REAL_QZ_H
12 
13 namespace Eigen {
14 
57  template<typename _MatrixType> class RealQZ
58  {
59  public:
60  typedef _MatrixType MatrixType;
61  enum {
62  RowsAtCompileTime = MatrixType::RowsAtCompileTime,
63  ColsAtCompileTime = MatrixType::ColsAtCompileTime,
64  Options = MatrixType::Options,
65  MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
66  MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67  };
68  typedef typename MatrixType::Scalar Scalar;
70  typedef Eigen::Index Index;
71 
74 
87  m_S(size, size),
88  m_T(size, size),
89  m_Q(size, size),
90  m_Z(size, size),
91  m_workspace(size*2),
92  m_maxIters(400),
93  m_isInitialized(false),
94  m_computeQZ(true)
95  {}
96 
105  RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true) :
106  m_S(A.rows(),A.cols()),
107  m_T(A.rows(),A.cols()),
108  m_Q(A.rows(),A.cols()),
109  m_Z(A.rows(),A.cols()),
110  m_workspace(A.rows()*2),
111  m_maxIters(400),
112  m_isInitialized(false),
113  m_computeQZ(true)
114  {
115  compute(A, B, computeQZ);
116  }
117 
122  const MatrixType& matrixQ() const {
123  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
124  eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
125  return m_Q;
126  }
127 
132  const MatrixType& matrixZ() const {
133  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
134  eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
135  return m_Z;
136  }
137 
142  const MatrixType& matrixS() const {
143  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
144  return m_S;
145  }
146 
151  const MatrixType& matrixT() const {
152  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
153  return m_T;
154  }
155 
163  RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);
164 
170  {
171  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
172  return m_info;
173  }
174 
177  Index iterations() const
178  {
179  eigen_assert(m_isInitialized && "RealQZ is not initialized.");
180  return m_global_iter;
181  }
182 
186  RealQZ& setMaxIterations(Index maxIters)
187  {
188  m_maxIters = maxIters;
189  return *this;
190  }
191 
192  private:
193 
194  MatrixType m_S, m_T, m_Q, m_Z;
197  Index m_maxIters;
202 
207 
208  void hessenbergTriangular();
209  void computeNorms();
210  Index findSmallSubdiagEntry(Index iu);
211  Index findSmallDiagEntry(Index f, Index l);
212  void splitOffTwoRows(Index i);
213  void pushDownZero(Index z, Index f, Index l);
214  void step(Index f, Index l, Index iter);
215 
216  }; // RealQZ
217 
219  template<typename MatrixType>
221  {
222 
223  const Index dim = m_S.cols();
224 
225  // perform QR decomposition of T, overwrite T with R, save Q
227  m_T = qrT.matrixQR();
228  m_T.template triangularView<StrictlyLower>().setZero();
229  m_Q = qrT.householderQ();
230  // overwrite S with Q* S
231  m_S.applyOnTheLeft(m_Q.adjoint());
232  // init Z as Identity
233  if (m_computeQZ)
234  m_Z = MatrixType::Identity(dim,dim);
235  // reduce S to upper Hessenberg with Givens rotations
236  for (Index j=0; j<=dim-3; j++) {
237  for (Index i=dim-1; i>=j+2; i--) {
238  JRs G;
239  // kill S(i,j)
240  if(m_S.coeff(i,j) != 0)
241  {
242  G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
243  m_S.coeffRef(i,j) = Scalar(0.0);
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());
246  // update Q
247  if (m_computeQZ)
248  m_Q.applyOnTheRight(i-1,i,G);
249  }
250  // kill T(i,i-1)
251  if(m_T.coeff(i,i-1)!=Scalar(0))
252  {
253  G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
254  m_T.coeffRef(i,i-1) = Scalar(0.0);
255  m_S.applyOnTheRight(i,i-1,G);
256  m_T.topRows(i).applyOnTheRight(i,i-1,G);
257  // update Z
258  if (m_computeQZ)
259  m_Z.applyOnTheLeft(i,i-1,G.adjoint());
260  }
261  }
262  }
263  }
264 
266  template<typename MatrixType>
268  {
269  const Index size = m_S.cols();
270  m_normOfS = Scalar(0.0);
271  m_normOfT = Scalar(0.0);
272  for (Index j = 0; j < size; ++j)
273  {
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();
276  }
277  }
278 
279 
281  template<typename MatrixType>
283  {
284  using std::abs;
285  Index res = iu;
286  while (res > 0)
287  {
288  Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
289  if (s == Scalar(0.0))
290  s = m_normOfS;
291  if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
292  break;
293  res--;
294  }
295  return res;
296  }
297 
299  template<typename MatrixType>
301  {
302  using std::abs;
303  Index res = l;
304  while (res >= f) {
305  if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
306  break;
307  res--;
308  }
309  return res;
310  }
311 
313  template<typename MatrixType>
315  {
316  using std::abs;
317  using std::sqrt;
318  const Index dim=m_S.cols();
319  if (abs(m_S.coeff(i+1,i))==Scalar(0))
320  return;
321  Index j = findSmallDiagEntry(i,i+1);
322  if (j==i-1)
323  {
324  // block of (S T^{-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));
327  Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1));
328  Scalar q = p*p + STi(1,0)*STi(0,1);
329  if (q>=0) {
330  Scalar z = sqrt(q);
331  // one QR-like iteration for ABi - lambda I
332  // is enough - when we know exact eigenvalue in advance,
333  // convergence is immediate
334  JRs G;
335  if (p>=0)
336  G.makeGivens(p + z, STi(1,0));
337  else
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());
341  // update Q
342  if (m_computeQZ)
343  m_Q.applyOnTheRight(i,i+1,G);
344 
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);
348  // update Z
349  if (m_computeQZ)
350  m_Z.applyOnTheLeft(i+1,i,G.adjoint());
351 
352  m_S.coeffRef(i+1,i) = Scalar(0.0);
353  m_T.coeffRef(i+1,i) = Scalar(0.0);
354  }
355  }
356  else
357  {
358  pushDownZero(j,i,i+1);
359  }
360  }
361 
363  template<typename MatrixType>
365  {
366  JRs G;
367  const Index dim = m_S.cols();
368  for (Index zz=z; zz<l; zz++)
369  {
370  // push 0 down
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);
376  // update Q
377  if (m_computeQZ)
378  m_Q.applyOnTheRight(zz,zz+1,G);
379  // kill S(zz+1, zz-1)
380  if (zz>f)
381  {
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);
386  // update Z
387  if (m_computeQZ)
388  m_Z.applyOnTheLeft(zz,zz-1,G.adjoint());
389  }
390  }
391  // finally kill S(l,l-1)
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);
395  m_S.coeffRef(l,l-1)=Scalar(0.0);
396  // update Z
397  if (m_computeQZ)
398  m_Z.applyOnTheLeft(l,l-1,G.adjoint());
399  }
400 
402  template<typename MatrixType>
404  {
405  using std::abs;
406  const Index dim = m_S.cols();
407 
408  // x, y, z
409  Scalar x, y, z;
410  if (iter==10)
411  {
412  // Wilkinson ad hoc shift
413  const Scalar
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);
423  Scalar ss = abs(a87*b77i) + abs(a98*b88i),
424  lpl = Scalar(1.5)*ss,
425  ll = ss*ss;
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;
431  }
432  else if (iter==16)
433  {
434  // another exceptional shift
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);
438  z = 0;
439  }
440  else if (iter>23 && !(iter%8))
441  {
442  // extremely exceptional shift
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);
446  }
447  else
448  {
449  // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1
450  // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where
451  // U and V are 2x2 bottom right sub matrices of A and B. Thus:
452  // = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1)
453  // = AB^-1AB^-1 + det(M) - tr(M)(AB^-1)
454  // Since we are only interested in having x, y, z with a correct ratio, we have:
455  const Scalar
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),
459 
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),
462 
463  b11 = m_T.coeff(f,f), b12 = m_T.coeff(f,f+1),
464  b22 = m_T.coeff(f+1,f+1),
465 
466  b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l),
467  b99 = m_T.coeff(l,l);
468 
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);
472  z = a32/b22;
473  }
474 
475  JRs G;
476 
477  for (Index k=f; k<=l-2; k++)
478  {
479  // variables for Householder reflections
480  Vector2s essential2;
481  Scalar tau, beta;
482 
483  Vector3s hr(x,y,z);
484 
485  // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1)
486  hr.makeHouseholderInPlace(tau, beta);
487  essential2 = hr.template bottomRows<2>();
488  Index fc=(std::max)(k-1,Index(0)); // first col to update
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());
491  if (m_computeQZ)
492  m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data());
493  if (k>f)
494  m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0);
495 
496  // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k)
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>();
500  {
501  Index lr = (std::min)(k+4,dim); // last row to update
503  // S
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();
508  // T
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();
513  }
514  if (m_computeQZ)
515  {
516  // Z
518  tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k));
519  tmp += m_Z.row(k+2);
520  m_Z.row(k+2) -= tau*tmp;
521  m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp);
522  }
523  m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0);
524 
525  // Z_{k2} to annihilate T(k+1,k)
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);
529  // update Z
530  if (m_computeQZ)
531  m_Z.applyOnTheLeft(k+1,k,G.adjoint());
532  m_T.coeffRef(k+1,k) = Scalar(0.0);
533 
534  // update x,y,z
535  x = m_S.coeff(k+1,k);
536  y = m_S.coeff(k+2,k);
537  if (k < l-2)
538  z = m_S.coeff(k+3,k);
539  } // loop over k
540 
541  // Q_{n-1} to annihilate y = S(l,l-2)
542  G.makeGivens(x,y);
543  m_S.applyOnTheLeft(l-1,l,G.adjoint());
544  m_T.applyOnTheLeft(l-1,l,G.adjoint());
545  if (m_computeQZ)
546  m_Q.applyOnTheRight(l-1,l,G);
547  m_S.coeffRef(l,l-2) = Scalar(0.0);
548 
549  // Z_{n-1} to annihilate T(l,l-1)
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);
553  if (m_computeQZ)
554  m_Z.applyOnTheLeft(l,l-1,G.adjoint());
555  m_T.coeffRef(l,l-1) = Scalar(0.0);
556  }
557 
558  template<typename MatrixType>
559  RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ)
560  {
561 
562  const Index dim = A_in.cols();
563 
564  eigen_assert (A_in.rows()==dim && A_in.cols()==dim
565  && B_in.rows()==dim && B_in.cols()==dim
566  && "Need square matrices of the same dimension");
567 
568  m_isInitialized = true;
569  m_computeQZ = computeQZ;
570  m_S = A_in; m_T = B_in;
571  m_workspace.resize(dim*2);
572  m_global_iter = 0;
573 
574  // entrance point: hessenberg triangular decomposition
576  // compute L1 vector norms of T, S into m_normOfS, m_normOfT
577  computeNorms();
578 
579  Index l = dim-1,
580  f,
581  local_iter = 0;
582 
583  while (l>0 && local_iter<m_maxIters)
584  {
585  f = findSmallSubdiagEntry(l);
586  // now rows and columns f..l (including) decouple from the rest of the problem
587  if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0);
588  if (f == l) // One root found
589  {
590  l--;
591  local_iter = 0;
592  }
593  else if (f == l-1) // Two roots found
594  {
595  splitOffTwoRows(f);
596  l -= 2;
597  local_iter = 0;
598  }
599  else // No convergence yet
600  {
601  // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations
602  Index z = findSmallDiagEntry(f,l);
603  if (z>=f)
604  {
605  // zero found
606  pushDownZero(z,f,l);
607  }
608  else
609  {
610  // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg
611  // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to
612  // apply a QR-like iteration to rows and columns f..l.
613  step(f,l, local_iter);
614  local_iter++;
615  m_global_iter++;
616  }
617  }
618  }
619  // check if we converged before reaching iterations limit
620  m_info = (local_iter<m_maxIters) ? Success : NoConvergence;
621 
622  // For each non triangular 2x2 diagonal block of S,
623  // reduce the respective 2x2 diagonal block of T to positive diagonal form using 2x2 SVD.
624  // This step is not mandatory for QZ, but it does help further extraction of eigenvalues/eigenvectors,
625  // and is in par with Lapack/Matlab QZ.
626  if(m_info==Success)
627  {
628  for(Index i=0; i<dim-1; ++i)
629  {
630  if(m_S.coeff(i+1, i) != Scalar(0))
631  {
632  JacobiRotation<Scalar> j_left, j_right;
633  internal::real_2x2_jacobi_svd(m_T, i, i+1, &j_left, &j_right);
634 
635  // Apply resulting Jacobi rotations
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);
640  m_T(i+1,i) = m_T(i,i+1) = Scalar(0);
641 
642  if(m_computeQZ) {
643  m_Q.applyOnTheRight(i,i+1,j_left.transpose());
644  m_Z.applyOnTheLeft(i,i+1,j_right.transpose());
645  }
646 
647  i++;
648  }
649  }
650  }
651 
652  return *this;
653  } // end compute
654 
655 } // end namespace Eigen
656 
657 #endif //EIGEN_REAL_QZ
Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > EigenvalueType
Definition: RealQZ.h:72
void pushDownZero(Index z, Index f, Index l)
Definition: RealQZ.h:364
Index m_maxIters
Definition: RealQZ.h:197
SCALAR Scalar
Definition: bench_gemm.cpp:46
#define max(a, b)
Definition: datatypes.h:20
Index m_global_iter
Definition: RealQZ.h:201
MatrixType m_Z
Definition: RealQZ.h:194
Scalar * y
const MatrixType & matrixS() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:142
std::complex< typename NumTraits< Scalar >::Real > ComplexScalar
Definition: RealQZ.h:69
MatrixType m_Q
Definition: RealQZ.h:194
JacobiRotation< float > G
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE const Scalar * data() const
Eigen::Index Index
Definition: RealQZ.h:70
Index iterations() const
Returns number of performed QR-like iterations.
Definition: RealQZ.h:177
#define min(a, b)
Definition: datatypes.h:19
Index findSmallSubdiagEntry(Index iu)
Definition: RealQZ.h:282
MatrixType m_T
Definition: RealQZ.h:194
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const MatrixType & matrixQR() const
RealQZ(Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
Default constructor.
Definition: RealQZ.h:86
Namespace containing all symbols from the Eigen library.
Definition: jet.h:637
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NRowsBlockXpr< internal::get_fixed_value< NRowsType >::value >::Type topRows(NRowsType n)
Definition: BlockMethods.h:570
Rotation given by a cosine-sine pair.
iterator iter(handle obj)
Definition: pytypes.h:2273
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:232
void step(Index f, Index l, Index iter)
Definition: RealQZ.h:403
JacobiRotation< Scalar > JRs
Definition: RealQZ.h:206
RealQZ & setMaxIterations(Index maxIters)
Definition: RealQZ.h:186
Matrix< Scalar, 3, 1 > Vector3s
Definition: RealQZ.h:203
cout<< "Here is the matrix m:"<< endl<< m<< endl;Matrix< ptrdiff_t, 3, 1 > res
Scalar Scalar int size
Definition: benchVecAdd.cpp:17
void real_2x2_jacobi_svd(const MatrixType &matrix, Index p, Index q, JacobiRotation< RealScalar > *j_left, JacobiRotation< RealScalar > *j_right)
Definition: RealSvd2x2.h:19
static const Line3 l(Rot3(), 1, 1)
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE void resize(Index rows, Index cols)
MatrixType m_S
Definition: RealQZ.h:194
bool m_computeQZ
Definition: RealQZ.h:199
EIGEN_DEVICE_FUNC void makeGivens(const Scalar &p, const Scalar &q, Scalar *r=0)
Definition: Jacobi.h:162
RealQZ(const MatrixType &A, const MatrixType &B, bool computeQZ=true)
Constructor; computes real QZ decomposition of given matrices.
Definition: RealQZ.h:105
MatrixType::Scalar Scalar
Definition: RealQZ.h:68
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
#define eigen_assert(x)
Definition: Macros.h:1037
Matrix< Scalar, Dynamic, 1 > m_workspace
Definition: RealQZ.h:195
Point2(* f)(const Point3 &, OptionalJacobian< 2, 3 >)
RealScalar s
EIGEN_DEVICE_FUNC const Scalar & q
bool m_isInitialized
Definition: RealQZ.h:198
HouseholderSequenceType householderQ() const
Matrix< Scalar, 2, 1 > Vector2s
Definition: RealQZ.h:204
const MatrixType & matrixZ() const
Returns matrix Z in the QZ decomposition.
Definition: RealQZ.h:132
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: RealQZ.h:169
Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 > ColumnVectorType
Definition: RealQZ.h:73
EIGEN_DEVICE_FUNC JacobiRotation adjoint() const
Definition: Jacobi.h:67
static std::stringstream ss
Definition: testBTree.cpp:31
Scalar m_normOfS
Definition: RealQZ.h:200
mp::number< mp::cpp_dec_float< 100 >, mp::et_on > Real
Matrix< Scalar, 2, 2 > Matrix2s
Definition: RealQZ.h:205
EIGEN_DEVICE_FUNC JacobiRotation transpose() const
Definition: Jacobi.h:63
const MatrixType & matrixT() const
Returns matrix S in the QZ decomposition.
Definition: RealQZ.h:151
void splitOffTwoRows(Index i)
Definition: RealQZ.h:314
float * p
RealQZ & compute(const MatrixType &A, const MatrixType &B, bool computeQZ=true)
Computes QZ decomposition of given matrix.
Definition: RealQZ.h:559
Index findSmallDiagEntry(Index f, Index l)
Definition: RealQZ.h:300
ComputationInfo m_info
Definition: RealQZ.h:196
Scalar m_normOfT
Definition: RealQZ.h:200
_MatrixType MatrixType
Definition: RealQZ.h:60
void hessenbergTriangular()
Definition: RealQZ.h:220
Jet< T, N > sqrt(const Jet< T, N > &f)
Definition: jet.h:418
const int Dynamic
Definition: Constants.h:22
set noclip points set clip one set noclip two set bar set border lt lw set xdata set ydata set zdata set x2data set y2data set boxwidth set dummy x
EIGEN_DEVICE_FUNC EIGEN_STRONG_INLINE NColsBlockXpr< internal::get_fixed_value< NColsType >::value >::Type rightCols(NColsType n)
Definition: BlockMethods.h:872
#define abs(x)
Definition: datatypes.h:17
ComputationInfo
Definition: Constants.h:440
const MatrixType & matrixQ() const
Returns matrix Q in the QZ decomposition.
Definition: RealQZ.h:122
std::ptrdiff_t j
void computeNorms()
Definition: RealQZ.h:267
Performs a real QZ decomposition of a pair of square matrices.
Definition: RealQZ.h:57
v setZero(3)


gtsam
Author(s):
autogenerated on Tue Jul 4 2023 02:35:30