// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2012 Alexey Korepanov <kaikaikai@yandex.ru>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.

#ifndef EIGEN_REAL_QZ_H
#define EIGEN_REAL_QZ_H

// IWYU pragma: private
#include "./InternalHeaderCheck.h"

namespace Eigen {

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 *
 * \class RealQZ
 *
 * \brief Performs a real QZ decomposition of a pair of square matrices
 *
 * \tparam 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: \f$ A = Q S Z \f$, \f$ B = Q T Z \f$ 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, \f$ U^{-1} = U^T \f$. 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 \f$ A - z B \f$ 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: \include RealQZ_compute.cpp
 * Output: \include RealQZ_compute.out
 *
 * \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.
 *
 * \sa class RealSchur, class ComplexSchur, class EigenSolver, class ComplexEigenSolver
 */

template <typename MatrixType_>
class RealQZ {
 public:
  typedef MatrixType_ MatrixType;
  enum {
    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
    Options = internal::traits<MatrixType>::Options,
    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
  };
  typedef typename MatrixType::Scalar Scalar;
  typedef std::complex<typename NumTraits<Scalar>::Real> ComplexScalar;
  typedef Eigen::Index Index;  ///< \deprecated since Eigen 3.3

  typedef Matrix<ComplexScalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> EigenvalueType;
  typedef Matrix<Scalar, ColsAtCompileTime, 1, Options & ~RowMajor, MaxColsAtCompileTime, 1> ColumnVectorType;

  /** \brief Default constructor.
   *
   * \param [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 \p size parameter is only
   * used as a hint. It is not an error to give a wrong \p size, but it may
   * impair performance.
   *
   * \sa compute() for an example.
   */
  explicit RealQZ(Index size = RowsAtCompileTime == Dynamic ? 1 : 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) {}

  /** \brief Constructor; computes real QZ decomposition of given matrices
   *
   * \param[in]  A          Matrix A.
   * \param[in]  B          Matrix B.
   * \param[in]  computeQZ  If false, A and Z are not computed.
   *
   * This constructor calls compute() to compute the QZ decomposition.
   */
  RealQZ(const MatrixType& A, const MatrixType& B, bool computeQZ = true)
      : 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);
  }

  /** \brief Returns matrix Q in the QZ decomposition.
   *
   * \returns A const reference to the matrix Q.
   */
  const MatrixType& matrixQ() const {
    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;
  }

  /** \brief Returns matrix Z in the QZ decomposition.
   *
   * \returns A const reference to the matrix Z.
   */
  const MatrixType& matrixZ() const {
    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;
  }

  /** \brief Returns matrix S in the QZ decomposition.
   *
   * \returns A const reference to the matrix S.
   */
  const MatrixType& matrixS() const {
    eigen_assert(m_isInitialized && "RealQZ is not initialized.");
    return m_S;
  }

  /** \brief Returns matrix S in the QZ decomposition.
   *
   * \returns A const reference to the matrix S.
   */
  const MatrixType& matrixT() const {
    eigen_assert(m_isInitialized && "RealQZ is not initialized.");
    return m_T;
  }

  /** \brief Computes QZ decomposition of given matrix.
   *
   * \param[in]  A          Matrix A.
   * \param[in]  B          Matrix B.
   * \param[in]  computeQZ  If false, A and Z are not computed.
   * \returns    Reference to \c *this
   */
  RealQZ& compute(const MatrixType& A, const MatrixType& B, bool computeQZ = true);

  /** \brief Reports whether previous computation was successful.
   *
   * \returns \c Success if computation was successful, \c NoConvergence otherwise.
   */
  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "RealQZ is not initialized.");
    return m_info;
  }

  /** \brief Returns number of performed QR-like iterations.
   */
  Index iterations() const {
    eigen_assert(m_isInitialized && "RealQZ is not initialized.");
    return m_global_iter;
  }

  /** Sets the maximal number of iterations allowed to converge to one eigenvalue
   * or decouple the problem.
   */
  RealQZ& setMaxIterations(Index maxIters) {
    m_maxIters = maxIters;
    return *this;
  }

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

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

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

};  // RealQZ

/** \internal Reduces S and T to upper Hessenberg - triangular form */
template <typename MatrixType>
void RealQZ<MatrixType>::hessenbergTriangular() {
  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());
      }
    }
  }
}

/** \internal Computes vector L1 norms of S and T when in Hessenberg-Triangular form already */
template <typename MatrixType>
inline void RealQZ<MatrixType>::computeNorms() {
  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();
  }
}

/** \internal Look for single small sub-diagonal element S(res, res-1) and return res (or 0) */
template <typename MatrixType>
inline Index RealQZ<MatrixType>::findSmallSubdiagEntry(Index iu) {
  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;
}

/** \internal Look for single small diagonal element T(res, res) for res between f and l, and return res (or f-1)  */
template <typename MatrixType>
inline Index RealQZ<MatrixType>::findSmallDiagEntry(Index f, Index l) {
  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;
}

/** \internal decouple 2x2 diagonal block in rows i, i+1 if eigenvalues are real */
template <typename MatrixType>
inline void RealQZ<MatrixType>::splitOffTwoRows(Index i) {
  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);
  }
}

/** \internal use zero in T(z,z) to zero S(l,l-1), working in block f..l */
template <typename MatrixType>
inline void RealQZ<MatrixType>::pushDownZero(Index z, Index f, Index l) {
  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());
}

/** \internal QR-like iterative step for block f..l */
template <typename MatrixType>
inline void RealQZ<MatrixType>::step(Index f, Index l, Index iter) {
  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);
}

template <typename MatrixType>
RealQZ<MatrixType>& RealQZ<MatrixType>::compute(const MatrixType& A_in, const MatrixType& B_in, bool computeQZ) {
  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

}  // end namespace Eigen

#endif  // EIGEN_REAL_QZ