Eigen::RealQZ< MatrixType_ > Class Template Reference

Performs a real QZ decomposition of a pair of square matrices. More...

Public Types
enum	{   RowsAtCompileTime ,   ColsAtCompileTime ,   Options ,   MaxRowsAtCompileTime ,   MaxColsAtCompileTime }
typedef Matrix< Scalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 >	ColumnVectorType
typedef std::complex< typename NumTraits< Scalar >::Real >	ComplexScalar
typedef Matrix< ComplexScalar, ColsAtCompileTime, 1, Options &~RowMajor, MaxColsAtCompileTime, 1 >	EigenvalueType
typedef Eigen::Index	Index
typedef MatrixType_	MatrixType
typedef MatrixType::Scalar	Scalar

Public Member Functions
RealQZ &	compute (const MatrixType &A, const MatrixType &B, bool computeQZ=true)
	Computes QZ decomposition of given matrix. More...
ComputationInfo	info () const
	Reports whether previous computation was successful. More...
Index	iterations () const
	Returns number of performed QR-like iterations. More...
const MatrixType &	matrixQ () const
	Returns matrix Q in the QZ decomposition. More...
const MatrixType &	matrixS () const
	Returns matrix S in the QZ decomposition. More...
const MatrixType &	matrixT () const
	Returns matrix S in the QZ decomposition. More...
const MatrixType &	matrixZ () const
	Returns matrix Z in the QZ decomposition. More...
	RealQZ (const MatrixType &A, const MatrixType &B, bool computeQZ=true)
	Constructor; computes real QZ decomposition of given matrices. More...
	RealQZ (Index size=RowsAtCompileTime==Dynamic ? 1 :RowsAtCompileTime)
	Default constructor. More...
RealQZ &	setMaxIterations (Index maxIters)

Private Types
typedef JacobiRotation< Scalar >	JRs
typedef Matrix< Scalar, 2, 2 >	Matrix2s
typedef Matrix< Scalar, 2, 1 >	Vector2s
typedef Matrix< Scalar, 3, 1 >	Vector3s

Private Member Functions
void	computeNorms ()
Index	findSmallDiagEntry (Index f, Index l)
Index	findSmallSubdiagEntry (Index iu)
void	hessenbergTriangular ()
void	pushDownZero (Index z, Index f, Index l)
void	splitOffTwoRows (Index i)
void	step (Index f, Index l, Index iter)

Private Attributes
bool	m_computeQZ
Index	m_global_iter
ComputationInfo	m_info
bool	m_isInitialized
Index	m_maxIters
Scalar	m_normOfS
Scalar	m_normOfT
MatrixType	m_Q
MatrixType	m_S
MatrixType	m_T
Matrix< Scalar, Dynamic, 1 >	m_workspace
MatrixType	m_Z

Detailed Description

template<typename MatrixType_>
class Eigen::RealQZ< MatrixType_ >

Performs a real QZ decomposition of a pair of square matrices.

This is defined in the Eigenvalues module.

#include <Eigen/Eigenvalues> 

Template Parameters

MatrixType_

the type of the matrix of which we are computing the real QZ decomposition; this is expected to be an instantiation of the Matrix class template.

Given a real square matrices A and B, this class computes the real QZ decomposition: \( A = Q S Z \), \( B = Q T Z \) where Q and Z are real orthogonal matrixes, T is upper-triangular matrix, and S is upper quasi-triangular matrix. An orthogonal matrix is a matrix whose inverse is equal to its transpose, \( U^{-1} = U^T \). A quasi-triangular matrix is a block-triangular matrix whose diagonal consists of 1-by-1 blocks and 2-by-2 blocks where further reduction is impossible due to complex eigenvalues.

The eigenvalues of the pencil \( A - z B \) can be obtained from 1x1 and 2x2 blocks on the diagonals of S and T.

Call the function compute() to compute the real QZ decomposition of a given pair of matrices. Alternatively, you can use the RealQZ(const MatrixType& B, const MatrixType& B, bool computeQZ) constructor which computes the real QZ decomposition at construction time. Once the decomposition is computed, you can use the matrixS(), matrixT(), matrixQ() and matrixZ() functions to retrieve the matrices S, T, Q and Z in the decomposition. If computeQZ==false, some time is saved by not computing matrices Q and Z.

Example:

MatrixXf A = MatrixXf::Random(4,4);
MatrixXf B = MatrixXf::Random(4,4);
RealQZ<MatrixXf> qz(4); // preallocate space for 4x4 matrices
qz.compute(A,B);  // A = Q S Z,  B = Q T Z
 
// print original matrices and result of decomposition
cout << "A:\n" << A << "\n" << "B:\n" << B << "\n";
cout << "S:\n" << qz.matrixS() << "\n" << "T:\n" << qz.matrixT() << "\n";
cout << "Q:\n" << qz.matrixQ() << "\n" << "Z:\n" << qz.matrixZ() << "\n";
 
// verify precision
cout << "\nErrors:"
  << "\n|A-QSZ|: " << (A-qz.matrixQ()*qz.matrixS()*qz.matrixZ()).norm()
  << ", |B-QTZ|: " << (B-qz.matrixQ()*qz.matrixT()*qz.matrixZ()).norm()
  << "\n|QQ* - I|: " << (qz.matrixQ()*qz.matrixQ().adjoint() - MatrixXf::Identity(4,4)).norm()
  << ", |ZZ* - I|: " << (qz.matrixZ()*qz.matrixZ().adjoint() - MatrixXf::Identity(4,4)).norm()
  << "\n";

Output:

A:
   0.68   0.823  -0.444   -0.27
 -0.211  -0.605   0.108  0.0268
  0.566   -0.33 -0.0452   0.904
  0.597   0.536   0.258   0.832
B:
 0.271 -0.967 -0.687  0.998
 0.435 -0.514 -0.198 -0.563
-0.717 -0.726  -0.74 0.0259
 0.214  0.608 -0.782  0.678
S:
-0.668   1.26 -0.598 0.0941
 0.317  -0.27 -0.279   0.64
     0      0 -0.398 -0.164
     0      0      0  -1.12
T:
 -1.55      0  0.342  -0.54
     0   1.01 -0.457  0.128
     0      0  -1.25  0.438
     0      0      0  0.746
Q:
-0.587 -0.138  0.552  0.576
  0.19 -0.208  0.761 -0.585
-0.292  0.918  0.152 -0.223
-0.731  -0.31 -0.306 -0.526
Z:
-0.0204   0.213   -0.78   0.588
 -0.961  -0.184   -0.14  -0.153
 -0.269   0.783   0.462    0.32
-0.0674  -0.555   0.398   0.727
 
Errors:
|A-QSZ|: 1.36e-06, |B-QTZ|: 1.83e-06
|QQ* - I|: 8.18e-07, |ZZ* - I|: 7.36e-07

Note

The implementation is based on the algorithm in "Matrix Computations" by Gene H. Golub and Charles F. Van Loan, and a paper "An algorithm for generalized eigenvalue problems" by C.B.Moler and G.W.Stewart.

See also

class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver

Definition at line 59 of file RealQZ.h.

Member Typedef Documentation

ColumnVectorType

template<typename MatrixType_ >

typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> Eigen::RealQZ< MatrixType_ >::ColumnVectorType

Definition at line 75 of file RealQZ.h.

ComplexScalar

template<typename MatrixType_ >

typedef std::complex<typename NumTraits<Scalar>::Real> Eigen::RealQZ< MatrixType_ >::ComplexScalar

Definition at line 71 of file RealQZ.h.

EigenvalueType

template<typename MatrixType_ >

typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> Eigen::RealQZ< MatrixType_ >::EigenvalueType

Definition at line 74 of file RealQZ.h.

Index

template<typename MatrixType_ >

typedef Eigen::Index Eigen::RealQZ< MatrixType_ >::Index

Deprecated:

since Eigen 3.3

Definition at line 72 of file RealQZ.h.

JRs

template<typename MatrixType_ >

typedef JacobiRotation<Scalar> Eigen::RealQZ< MatrixType_ >::JRs

private

Definition at line 208 of file RealQZ.h.

Matrix2s

template<typename MatrixType_ >

typedef Matrix<Scalar,2,2> Eigen::RealQZ< MatrixType_ >::Matrix2s

private

Definition at line 207 of file RealQZ.h.

MatrixType

template<typename MatrixType_ >

typedef MatrixType_ Eigen::RealQZ< MatrixType_ >::MatrixType

Definition at line 62 of file RealQZ.h.

Scalar

template<typename MatrixType_ >

typedef MatrixType::Scalar Eigen::RealQZ< MatrixType_ >::Scalar

Definition at line 70 of file RealQZ.h.

Vector2s

template<typename MatrixType_ >

typedef Matrix<Scalar,2,1> Eigen::RealQZ< MatrixType_ >::Vector2s

private

Definition at line 206 of file RealQZ.h.

Vector3s

template<typename MatrixType_ >

typedef Matrix<Scalar,3,1> Eigen::RealQZ< MatrixType_ >::Vector3s

private

Definition at line 205 of file RealQZ.h.

Member Enumeration Documentation

anonymous enum

template<typename MatrixType_ >

anonymous enum

Enumerator
RowsAtCompileTime
ColsAtCompileTime
Options
MaxRowsAtCompileTime
MaxColsAtCompileTime

Definition at line 63 of file RealQZ.h.

           {
        RowsAtCompileTime = MatrixType::RowsAtCompileTime,
        ColsAtCompileTime = MatrixType::ColsAtCompileTime,
        Options = MatrixType::Options,
        MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
        MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
      };

Constructor & Destructor Documentation

RealQZ() [1/2]

template<typename MatrixType_ >

Eigen::RealQZ< MatrixType_ >::RealQZ ( Index  size = RowsAtCompileTime==Dynamic ? 1 : RowsAtCompileTime )

inlineexplicit

Default constructor.

Parameters

[in]

size

Positive integer, size of the matrix whose QZ decomposition will be computed.

The default constructor is useful in cases in which the user intends to perform decompositions via compute(). The size parameter is only used as a hint. It is not an error to give a wrong size, but it may impair performance.

See also

compute() for an example.

Definition at line 88 of file RealQZ.h.

                                                                  : RowsAtCompileTime) :
        m_S(size, size),
        m_T(size, size),
        m_Q(size, size),
        m_Z(size, size),
        m_workspace(size*2),
        m_maxIters(400),
        m_isInitialized(false),
        m_computeQZ(true)
      {}

RealQZ() [2/2]

template<typename MatrixType_ >

Eigen::RealQZ< MatrixType_ >::RealQZ ( const MatrixType &  A, const MatrixType &  B, bool  computeQZ = true  )

inline

Constructor; computes real QZ decomposition of given matrices.

Parameters

[in]	A	Matrix A.
[in]	B	Matrix B.
[in]	computeQZ	If false, A and Z are not computed.

This constructor calls compute() to compute the QZ decomposition.

Definition at line 107 of file RealQZ.h.

                                                                              :
        m_S(A.rows(),A.cols()),
        m_T(A.rows(),A.cols()),
        m_Q(A.rows(),A.cols()),
        m_Z(A.rows(),A.cols()),
        m_workspace(A.rows()*2),
        m_maxIters(400),
        m_isInitialized(false),
        m_computeQZ(true)
      {
        compute(A, B, computeQZ);
      }

Member Function Documentation

compute()

template<typename MatrixType >

RealQZ< MatrixType > & Eigen::RealQZ< MatrixType >::compute	(	const MatrixType &	A,
		const MatrixType &	B,
		bool	computeQZ = true
	)

Computes QZ decomposition of given matrix.

Parameters

[in]	A	Matrix A.
[in]	B	Matrix B.
[in]	computeQZ	If false, A and Z are not computed.

Returns

Reference to *this

Definition at line 561 of file RealQZ.h.

    {
 
      const Index dim = A_in.cols();
 
      eigen_assert (A_in.rows()==dim && A_in.cols()==dim 
          && B_in.rows()==dim && B_in.cols()==dim 
          && "Need square matrices of the same dimension");
 
      m_isInitialized = true;
      m_computeQZ = computeQZ;
      m_S = A_in; m_T = B_in;
      m_workspace.resize(dim*2);
      m_global_iter = 0;
 
      // entrance point: hessenberg triangular decomposition
      hessenbergTriangular();
      // compute L1 vector norms of T, S into m_normOfS, m_normOfT
      computeNorms();
 
      Index l = dim-1, 
            f, 
            local_iter = 0;
 
      while (l>0 && local_iter<m_maxIters)
      {
        f = findSmallSubdiagEntry(l);
        // now rows and columns f..l (including) decouple from the rest of the problem
        if (f>0) m_S.coeffRef(f,f-1) = Scalar(0.0);
        if (f == l) // One root found
        {
          l--;
          local_iter = 0;
        }
        else if (f == l-1) // Two roots found
        {
          splitOffTwoRows(f);
          l -= 2;
          local_iter = 0;
        }
        else // No convergence yet
        {
          // if there's zero on diagonal of T, we can isolate an eigenvalue with Givens rotations
          Index z = findSmallDiagEntry(f,l);
          if (z>=f)
          {
            // zero found
            pushDownZero(z,f,l);
          }
          else
          {
            // We are sure now that S.block(f,f, l-f+1,l-f+1) is underuced upper-Hessenberg 
            // and T.block(f,f, l-f+1,l-f+1) is invertible uper-triangular, which allows to
            // apply a QR-like iteration to rows and columns f..l.
            step(f,l, local_iter);
            local_iter++;
            m_global_iter++;
          }
        }
      }
      // 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(!numext::is_exactly_zero(m_S.coeff(i + 1, i)))
          {
            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_S.applyOnTheLeft(i,i+1,j_left);
            m_S.applyOnTheRight(i,i+1,j_right);
            m_T.applyOnTheLeft(i,i+1,j_left);
            m_T.applyOnTheRight(i,i+1,j_right);
            m_T(i+1,i) = 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

computeNorms()

template<typename MatrixType >

void Eigen::RealQZ< MatrixType >::computeNorms

inlineprivate

Definition at line 269 of file RealQZ.h.

    {
      const Index size = m_S.cols();
      m_normOfS = Scalar(0.0);
      m_normOfT = Scalar(0.0);
      for (Index j = 0; j < size; ++j)
      {
        m_normOfS += m_S.col(j).segment(0, (std::min)(size,j+2)).cwiseAbs().sum();
        m_normOfT += m_T.row(j).segment(j, size - j).cwiseAbs().sum();
      }
    }

findSmallDiagEntry()

template<typename MatrixType >

Index Eigen::RealQZ< MatrixType >::findSmallDiagEntry ( Index  f, Index  l  )

inlineprivate

Definition at line 302 of file RealQZ.h.

    {
      using std::abs;
      Index res = l;
      while (res >= f) {
        if (abs(m_T.coeff(res,res)) <= NumTraits<Scalar>::epsilon() * m_normOfT)
          break;
        res--;
      }
      return res;
    }

findSmallSubdiagEntry()

template<typename MatrixType >

Index Eigen::RealQZ< MatrixType >::findSmallSubdiagEntry ( Index  iu )

inlineprivate

Definition at line 284 of file RealQZ.h.

    {
      using std::abs;
      Index res = iu;
      while (res > 0)
      {
        Scalar s = abs(m_S.coeff(res-1,res-1)) + abs(m_S.coeff(res,res));
        if (numext::is_exactly_zero(s))
          s = m_normOfS;
        if (abs(m_S.coeff(res,res-1)) < NumTraits<Scalar>::epsilon() * s)
          break;
        res--;
      }
      return res;
    }

hessenbergTriangular()

template<typename MatrixType >

void Eigen::RealQZ< MatrixType >::hessenbergTriangular

private

Definition at line 222 of file RealQZ.h.

    {
 
      const Index dim = m_S.cols();
 
      // perform QR decomposition of T, overwrite T with R, save Q
      HouseholderQR<MatrixType> qrT(m_T);
      m_T = qrT.matrixQR();
      m_T.template triangularView<StrictlyLower>().setZero();
      m_Q = qrT.householderQ();
      // overwrite S with Q* S
      m_S.applyOnTheLeft(m_Q.adjoint());
      // init Z as Identity
      if (m_computeQZ)
        m_Z = MatrixType::Identity(dim,dim);
      // reduce S to upper Hessenberg with Givens rotations
      for (Index j=0; j<=dim-3; j++) {
        for (Index i=dim-1; i>=j+2; i--) {
          JRs G;
          // kill S(i,j)
          if(!numext::is_exactly_zero(m_S.coeff(i, j)))
          {
            G.makeGivens(m_S.coeff(i-1,j), m_S.coeff(i,j), &m_S.coeffRef(i-1, j));
            m_S.coeffRef(i,j) = Scalar(0.0);
            m_S.rightCols(dim-j-1).applyOnTheLeft(i-1,i,G.adjoint());
            m_T.rightCols(dim-i+1).applyOnTheLeft(i-1,i,G.adjoint());
            // update Q
            if (m_computeQZ)
              m_Q.applyOnTheRight(i-1,i,G);
          }
          // kill T(i,i-1)
          if(!numext::is_exactly_zero(m_T.coeff(i, i - 1)))
          {
            G.makeGivens(m_T.coeff(i,i), m_T.coeff(i,i-1), &m_T.coeffRef(i,i));
            m_T.coeffRef(i,i-1) = Scalar(0.0);
            m_S.applyOnTheRight(i,i-1,G);
            m_T.topRows(i).applyOnTheRight(i,i-1,G);
            // update Z
            if (m_computeQZ)
              m_Z.applyOnTheLeft(i,i-1,G.adjoint());
          }
        }
      }
    }

info()

template<typename MatrixType_ >

ComputationInfo Eigen::RealQZ< MatrixType_ >::info ( ) const

inline

Reports whether previous computation was successful.

Returns

Success if computation was successful, NoConvergence otherwise.

Definition at line 171 of file RealQZ.h.

      {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        return m_info;
      }

iterations()

template<typename MatrixType_ >

Index Eigen::RealQZ< MatrixType_ >::iterations ( ) const

inline

Returns number of performed QR-like iterations.

Definition at line 179 of file RealQZ.h.

      {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        return m_global_iter;
      }

matrixQ()

template<typename MatrixType_ >

const MatrixType& Eigen::RealQZ< MatrixType_ >::matrixQ ( ) const

inline

Returns matrix Q in the QZ decomposition.

Returns

A const reference to the matrix Q.

Definition at line 124 of file RealQZ.h.

                                        {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
        return m_Q;
      }

matrixS()

template<typename MatrixType_ >

const MatrixType& Eigen::RealQZ< MatrixType_ >::matrixS ( ) const

inline

Returns matrix S in the QZ decomposition.

Returns

A const reference to the matrix S.

Definition at line 144 of file RealQZ.h.

                                        {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        return m_S;
      }

matrixT()

template<typename MatrixType_ >

const MatrixType& Eigen::RealQZ< MatrixType_ >::matrixT ( ) const

inline

Returns matrix S in the QZ decomposition.

Returns

A const reference to the matrix S.

Definition at line 153 of file RealQZ.h.

                                        {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        return m_T;
      }

matrixZ()

template<typename MatrixType_ >

const MatrixType& Eigen::RealQZ< MatrixType_ >::matrixZ ( ) const

inline

Returns matrix Z in the QZ decomposition.

Returns

A const reference to the matrix Z.

Definition at line 134 of file RealQZ.h.

                                        {
        eigen_assert(m_isInitialized && "RealQZ is not initialized.");
        eigen_assert(m_computeQZ && "The matrices Q and Z have not been computed during the QZ decomposition.");
        return m_Z;
      }

pushDownZero()

template<typename MatrixType >

void Eigen::RealQZ< MatrixType >::pushDownZero ( Index  z, Index  f, Index  l  )

inlineprivate

Definition at line 366 of file RealQZ.h.

    {
      JRs G;
      const Index dim = m_S.cols();
      for (Index zz=z; zz<l; zz++)
      {
        // push 0 down
        Index firstColS = zz>f ? (zz-1) : zz;
        G.makeGivens(m_T.coeff(zz, zz+1), m_T.coeff(zz+1, zz+1));
        m_S.rightCols(dim-firstColS).applyOnTheLeft(zz,zz+1,G.adjoint());
        m_T.rightCols(dim-zz).applyOnTheLeft(zz,zz+1,G.adjoint());
        m_T.coeffRef(zz+1,zz+1) = Scalar(0.0);
        // update Q
        if (m_computeQZ)
          m_Q.applyOnTheRight(zz,zz+1,G);
        // kill S(zz+1, zz-1)
        if (zz>f)
        {
          G.makeGivens(m_S.coeff(zz+1, zz), m_S.coeff(zz+1,zz-1));
          m_S.topRows(zz+2).applyOnTheRight(zz, zz-1,G);
          m_T.topRows(zz+1).applyOnTheRight(zz, zz-1,G);
          m_S.coeffRef(zz+1,zz-1) = Scalar(0.0);
          // update Z
          if (m_computeQZ)
            m_Z.applyOnTheLeft(zz,zz-1,G.adjoint());
        }
      }
      // finally kill S(l,l-1)
      G.makeGivens(m_S.coeff(l,l), m_S.coeff(l,l-1));
      m_S.applyOnTheRight(l,l-1,G);
      m_T.applyOnTheRight(l,l-1,G);
      m_S.coeffRef(l,l-1)=Scalar(0.0);
      // update Z
      if (m_computeQZ)
        m_Z.applyOnTheLeft(l,l-1,G.adjoint());
    }

setMaxIterations()

template<typename MatrixType_ >

RealQZ& Eigen::RealQZ< MatrixType_ >::setMaxIterations ( Index  maxIters )

inline

Sets the maximal number of iterations allowed to converge to one eigenvalue or decouple the problem.

Definition at line 188 of file RealQZ.h.

      {
        m_maxIters = maxIters;
        return *this;
      }

splitOffTwoRows()

template<typename MatrixType >

void Eigen::RealQZ< MatrixType >::splitOffTwoRows ( Index  i )

inlineprivate

Definition at line 316 of file RealQZ.h.

    {
      using std::abs;
      using std::sqrt;
      const Index dim=m_S.cols();
      if (numext::is_exactly_zero(abs(m_S.coeff(i + 1, i))))
        return;
      Index j = findSmallDiagEntry(i,i+1);
      if (j==i-1)
      {
        // block of (S T^{-1})
        Matrix2s STi = m_T.template block<2,2>(i,i).template triangularView<Upper>().
          template solve<OnTheRight>(m_S.template block<2,2>(i,i));
        Scalar p = Scalar(0.5)*(STi(0,0)-STi(1,1));
        Scalar q = p*p + STi(1,0)*STi(0,1);
        if (q>=0) {
          Scalar z = sqrt(q);
          // one QR-like iteration for ABi - lambda I
          // is enough - when we know exact eigenvalue in advance,
          // convergence is immediate
          JRs G;
          if (p>=0)
            G.makeGivens(p + z, STi(1,0));
          else
            G.makeGivens(p - z, STi(1,0));
          m_S.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
          m_T.rightCols(dim-i).applyOnTheLeft(i,i+1,G.adjoint());
          // update Q
          if (m_computeQZ)
            m_Q.applyOnTheRight(i,i+1,G);
 
          G.makeGivens(m_T.coeff(i+1,i+1), m_T.coeff(i+1,i));
          m_S.topRows(i+2).applyOnTheRight(i+1,i,G);
          m_T.topRows(i+2).applyOnTheRight(i+1,i,G);
          // update Z
          if (m_computeQZ)
            m_Z.applyOnTheLeft(i+1,i,G.adjoint());
 
          m_S.coeffRef(i+1,i) = Scalar(0.0);
          m_T.coeffRef(i+1,i) = Scalar(0.0);
        }
      }
      else
      {
        pushDownZero(j,i,i+1);
      }
    }

step()

template<typename MatrixType >

void Eigen::RealQZ< MatrixType >::step ( Index  f, Index  l, Index  iter  )

inlineprivate

Definition at line 405 of file RealQZ.h.

    {
      using std::abs;
      const Index dim = m_S.cols();
 
      // x, y, z
      Scalar x, y, z;
      if (iter==10)
      {
        // Wilkinson ad hoc shift
        const Scalar
          a11=m_S.coeff(f+0,f+0), a12=m_S.coeff(f+0,f+1),
          a21=m_S.coeff(f+1,f+0), a22=m_S.coeff(f+1,f+1), a32=m_S.coeff(f+2,f+1),
          b12=m_T.coeff(f+0,f+1),
          b11i=Scalar(1.0)/m_T.coeff(f+0,f+0),
          b22i=Scalar(1.0)/m_T.coeff(f+1,f+1),
          a87=m_S.coeff(l-1,l-2),
          a98=m_S.coeff(l-0,l-1),
          b77i=Scalar(1.0)/m_T.coeff(l-2,l-2),
          b88i=Scalar(1.0)/m_T.coeff(l-1,l-1);
        Scalar ss = abs(a87*b77i) + abs(a98*b88i),
               lpl = Scalar(1.5)*ss,
               ll = ss*ss;
        x = ll + a11*a11*b11i*b11i - lpl*a11*b11i + a12*a21*b11i*b22i
          - a11*a21*b12*b11i*b11i*b22i;
        y = a11*a21*b11i*b11i - lpl*a21*b11i + a21*a22*b11i*b22i 
          - a21*a21*b12*b11i*b11i*b22i;
        z = a21*a32*b11i*b22i;
      }
      else if (iter==16)
      {
        // another exceptional shift
        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) /
          (m_T.coeff(l-1,l-1)*m_T.coeff(l,l));
        y = m_S.coeff(f+1,f)/m_T.coeff(f,f);
        z = 0;
      }
      else if (iter>23 && !(iter%8))
      {
        // extremely exceptional shift
        x = internal::random<Scalar>(-1.0,1.0);
        y = internal::random<Scalar>(-1.0,1.0);
        z = internal::random<Scalar>(-1.0,1.0);
      }
      else
      {
        // Compute the shifts: (x,y,z,0...) = (AB^-1 - l1 I) (AB^-1 - l2 I) e1
        // where l1 and l2 are the eigenvalues of the 2x2 matrix C = U V^-1 where
        // U and V are 2x2 bottom right sub matrices of A and B. Thus:
        //  = AB^-1AB^-1 + l1 l2 I - (l1+l2)(AB^-1)
        //  = AB^-1AB^-1 + det(M) - tr(M)(AB^-1)
        // Since we are only interested in having x, y, z with a correct ratio, we have:
        const Scalar
          a11 = m_S.coeff(f,f),     a12 = m_S.coeff(f,f+1),
          a21 = m_S.coeff(f+1,f),   a22 = m_S.coeff(f+1,f+1),
                                    a32 = m_S.coeff(f+2,f+1),
 
          a88 = m_S.coeff(l-1,l-1), a89 = m_S.coeff(l-1,l),
          a98 = m_S.coeff(l,l-1),   a99 = m_S.coeff(l,l),
 
          b11 = m_T.coeff(f,f),     b12 = m_T.coeff(f,f+1),
                                    b22 = m_T.coeff(f+1,f+1),
 
          b88 = m_T.coeff(l-1,l-1), b89 = m_T.coeff(l-1,l),
                                    b99 = m_T.coeff(l,l);
 
        x = ( (a88/b88 - a11/b11)*(a99/b99 - a11/b11) - (a89/b99)*(a98/b88) + (a98/b88)*(b89/b99)*(a11/b11) ) * (b11/a21)
          + a12/b22 - (a11/b11)*(b12/b22);
        y = (a22/b22-a11/b11) - (a21/b11)*(b12/b22) - (a88/b88-a11/b11) - (a99/b99-a11/b11) + (a98/b88)*(b89/b99);
        z = a32/b22;
      }
 
      JRs G;
 
      for (Index k=f; k<=l-2; k++)
      {
        // variables for Householder reflections
        Vector2s essential2;
        Scalar tau, beta;
 
        Vector3s hr(x,y,z);
 
        // Q_k to annihilate S(k+1,k-1) and S(k+2,k-1)
        hr.makeHouseholderInPlace(tau, beta);
        essential2 = hr.template bottomRows<2>();
        Index fc=(std::max)(k-1,Index(0));  // first col to update
        m_S.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
        m_T.template middleRows<3>(k).rightCols(dim-fc).applyHouseholderOnTheLeft(essential2, tau, m_workspace.data());
        if (m_computeQZ)
          m_Q.template middleCols<3>(k).applyHouseholderOnTheRight(essential2, tau, m_workspace.data());
        if (k>f)
          m_S.coeffRef(k+2,k-1) = m_S.coeffRef(k+1,k-1) = Scalar(0.0);
 
        // Z_{k1} to annihilate T(k+2,k+1) and T(k+2,k)
        hr << m_T.coeff(k+2,k+2),m_T.coeff(k+2,k),m_T.coeff(k+2,k+1);
        hr.makeHouseholderInPlace(tau, beta);
        essential2 = hr.template bottomRows<2>();
        {
          Index lr = (std::min)(k+4,dim); // last row to update
          Map<Matrix<Scalar,Dynamic,1> > tmp(m_workspace.data(),lr);
          // S
          tmp = m_S.template middleCols<2>(k).topRows(lr) * essential2;
          tmp += m_S.col(k+2).head(lr);
          m_S.col(k+2).head(lr) -= tau*tmp;
          m_S.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
          // T
          tmp = m_T.template middleCols<2>(k).topRows(lr) * essential2;
          tmp += m_T.col(k+2).head(lr);
          m_T.col(k+2).head(lr) -= tau*tmp;
          m_T.template middleCols<2>(k).topRows(lr) -= (tau*tmp) * essential2.adjoint();
        }
        if (m_computeQZ)
        {
          // Z
          Map<Matrix<Scalar,1,Dynamic> > tmp(m_workspace.data(),dim);
          tmp = essential2.adjoint()*(m_Z.template middleRows<2>(k));
          tmp += m_Z.row(k+2);
          m_Z.row(k+2) -= tau*tmp;
          m_Z.template middleRows<2>(k) -= essential2 * (tau*tmp);
        }
        m_T.coeffRef(k+2,k) = m_T.coeffRef(k+2,k+1) = Scalar(0.0);
 
        // Z_{k2} to annihilate T(k+1,k)
        G.makeGivens(m_T.coeff(k+1,k+1), m_T.coeff(k+1,k));
        m_S.applyOnTheRight(k+1,k,G);
        m_T.applyOnTheRight(k+1,k,G);
        // update Z
        if (m_computeQZ)
          m_Z.applyOnTheLeft(k+1,k,G.adjoint());
        m_T.coeffRef(k+1,k) = Scalar(0.0);
 
        // update x,y,z
        x = m_S.coeff(k+1,k);
        y = m_S.coeff(k+2,k);
        if (k < l-2)
          z = m_S.coeff(k+3,k);
      } // loop over k
 
      // Q_{n-1} to annihilate y = S(l,l-2)
      G.makeGivens(x,y);
      m_S.applyOnTheLeft(l-1,l,G.adjoint());
      m_T.applyOnTheLeft(l-1,l,G.adjoint());
      if (m_computeQZ)
        m_Q.applyOnTheRight(l-1,l,G);
      m_S.coeffRef(l,l-2) = Scalar(0.0);
 
      // Z_{n-1} to annihilate T(l,l-1)
      G.makeGivens(m_T.coeff(l,l),m_T.coeff(l,l-1));
      m_S.applyOnTheRight(l,l-1,G);
      m_T.applyOnTheRight(l,l-1,G);
      if (m_computeQZ)
        m_Z.applyOnTheLeft(l,l-1,G.adjoint());
      m_T.coeffRef(l,l-1) = Scalar(0.0);
    }

Member Data Documentation

m_computeQZ

template<typename MatrixType_ >

bool Eigen::RealQZ< MatrixType_ >::m_computeQZ

private

Definition at line 201 of file RealQZ.h.

m_global_iter

template<typename MatrixType_ >

Index Eigen::RealQZ< MatrixType_ >::m_global_iter

private

Definition at line 203 of file RealQZ.h.

m_info

template<typename MatrixType_ >

ComputationInfo Eigen::RealQZ< MatrixType_ >::m_info

private

Definition at line 198 of file RealQZ.h.

m_isInitialized

template<typename MatrixType_ >

bool Eigen::RealQZ< MatrixType_ >::m_isInitialized

private

Definition at line 200 of file RealQZ.h.

m_maxIters

template<typename MatrixType_ >

Index Eigen::RealQZ< MatrixType_ >::m_maxIters

private

Definition at line 199 of file RealQZ.h.

m_normOfS

template<typename MatrixType_ >

Scalar Eigen::RealQZ< MatrixType_ >::m_normOfS

private

Definition at line 202 of file RealQZ.h.

m_normOfT

template<typename MatrixType_ >

Scalar Eigen::RealQZ< MatrixType_ >::m_normOfT

private

Definition at line 202 of file RealQZ.h.

m_Q

template<typename MatrixType_ >

MatrixType Eigen::RealQZ< MatrixType_ >::m_Q

private

Definition at line 196 of file RealQZ.h.

m_S

template<typename MatrixType_ >

MatrixType Eigen::RealQZ< MatrixType_ >::m_S

private

Definition at line 196 of file RealQZ.h.

m_T

template<typename MatrixType_ >

MatrixType Eigen::RealQZ< MatrixType_ >::m_T

private

Definition at line 196 of file RealQZ.h.

m_workspace

template<typename MatrixType_ >

Matrix<Scalar,Dynamic,1> Eigen::RealQZ< MatrixType_ >::m_workspace

private

Definition at line 197 of file RealQZ.h.

m_Z

template<typename MatrixType_ >

MatrixType Eigen::RealQZ< MatrixType_ >::m_Z

private

Definition at line 196 of file RealQZ.h.

---

The documentation for this class was generated from the following file:

• RealQZ.h