ComplexSchur.h

// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2009 Claire Maurice
// Copyright (C) 2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010,2012 Jitse Niesen <jitse@maths.leeds.ac.uk>
//
// 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/.
// SPDX-License-Identifier: MPL-2.0
 
#ifndef EIGEN_COMPLEX_SCHUR_H
#define EIGEN_COMPLEX_SCHUR_H
 
#include "./HessenbergDecomposition.h"
 
// IWYU pragma: private
#include "./InternalHeaderCheck.h"
 
namespace Eigen {
 
namespace internal {
template <typename MatrixType, bool IsComplex>
struct complex_schur_reduce_to_hessenberg;
}
 
template <typename MatrixType_>
class ComplexEigenSolver;

template <typename MatrixType_>
class ComplexSchur {
 public:
  using MatrixType = MatrixType_;
  enum {
    RowsAtCompileTime = MatrixType::RowsAtCompileTime,
    ColsAtCompileTime = MatrixType::ColsAtCompileTime,
    Options = internal::plain_object_options<MatrixType>::value,
    MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
    MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
  };

  using Scalar = typename MatrixType::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  using Index = Eigen::Index;  

  using ComplexScalar = internal::make_complex_t<Scalar>;

  using ComplexMatrixType =
      Matrix<ComplexScalar, RowsAtCompileTime, ColsAtCompileTime, Options, MaxRowsAtCompileTime, MaxColsAtCompileTime>;

  using MatrixTType = std::conditional_t<internal::is_ref<MatrixType>::value, MatrixType, ComplexMatrixType>;

  explicit ComplexSchur(Index size = RowsAtCompileTime == Dynamic ? 1 : RowsAtCompileTime)
      : m_matT(size, size),
        m_matU(size, size),
        m_hess(size),
        m_isInitialized(false),
        m_matUisUptodate(false),
        m_maxIters(-1) {}

  template <typename InputType>
  explicit ComplexSchur(const EigenBase<InputType>& matrix, bool computeU = true)
      : m_matT(matrix.rows(), matrix.cols()),
        m_matU(matrix.rows(), matrix.cols()),
        m_hess(matrix.rows()),
        m_isInitialized(false),
        m_matUisUptodate(false),
        m_maxIters(-1) {
    compute(matrix.derived(), computeU);
  }

  template <typename InputType, bool IsRef = internal::is_ref<MatrixType>::value, std::enable_if_t<IsRef, int> = 0>
  explicit ComplexSchur(EigenBase<InputType>& matrix, bool computeU = true)
      : m_matT(matrix.derived()),
        m_matU(matrix.rows(), matrix.cols()),
        m_hess(matrix.rows()),
        m_isInitialized(false),
        m_matUisUptodate(false),
        m_maxIters(-1) {
    computeInPlace(computeU);
  }

  const ComplexMatrixType& matrixU() const {
    eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
    eigen_assert(m_matUisUptodate && "The matrix U has not been computed during the ComplexSchur decomposition.");
    return m_matU;
  }

  const MatrixTType& matrixT() const {
    eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
    return m_matT;
  }

  template <typename InputType>
  ComplexSchur& compute(const EigenBase<InputType>& matrix, bool computeU = true);

  template <typename HessMatrixType, typename OrthMatrixType>
  ComplexSchur& computeFromHessenberg(const HessMatrixType& matrixH, const OrthMatrixType& matrixQ,
                                      bool computeU = true);

  ComputationInfo info() const {
    eigen_assert(m_isInitialized && "ComplexSchur is not initialized.");
    return m_info;
  }

  ComplexSchur& setMaxIterations(Index maxIters) {
    m_maxIters = maxIters;
    return *this;
  }

  Index getMaxIterations() const { return m_maxIters; }

  static const int m_maxIterationsPerRow = 30;
 
 protected:
  EIGEN_STATIC_ASSERT(NumTraits<Scalar>::IsComplex || !internal::is_ref<MatrixType>::value,
                      INPLACE_COMPLEXSCHUR_REQUIRES_A_COMPLEX_MATRIX_TYPE)
 
  MatrixTType m_matT;
  ComplexMatrixType m_matU;
  internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex> m_hess;
  ComputationInfo m_info;
  bool m_isInitialized;
  bool m_matUisUptodate;
  Index m_maxIters;
 
 private:
  // ComplexEigenSolver computes the eigenvectors of T within this storage.
  friend class ComplexEigenSolver<MatrixType>;
 
  ComplexSchur& computeInPlace(bool computeU);
  bool subdiagonalEntryIsNeglegible(Index i);
  ComplexScalar computeShift(Index iu, Index iter);
  void reduceToTriangularForm(bool computeU);
  friend struct internal::complex_schur_reduce_to_hessenberg<MatrixType, NumTraits<Scalar>::IsComplex>;
};

template <typename MatrixType>
inline bool ComplexSchur<MatrixType>::subdiagonalEntryIsNeglegible(Index i) {
  RealScalar d = numext::norm1(m_matT.coeff(i, i)) + numext::norm1(m_matT.coeff(i + 1, i + 1));
  RealScalar sd = numext::norm1(m_matT.coeff(i + 1, i));
  if (internal::isMuchSmallerThan(sd, d, NumTraits<RealScalar>::epsilon())) {
    m_matT.coeffRef(i + 1, i) = ComplexScalar(0);
    return true;
  }
  return false;
}

template <typename MatrixType>
typename ComplexSchur<MatrixType>::ComplexScalar ComplexSchur<MatrixType>::computeShift(Index iu, Index iter) {
  using std::abs;
  if ((iter == 10 || iter == 20) && iu > 1) {
    // exceptional shift, taken from http://www.netlib.org/eispack/comqr.f
    return ComplexSchur<MatrixType>::ComplexScalar(abs(numext::real(m_matT.coeff(iu, iu - 1))) +
                                                   abs(numext::real(m_matT.coeff(iu - 1, iu - 2))));
  }
 
  // compute the shift as one of the eigenvalues of t, the 2x2
  // diagonal block on the bottom of the active submatrix
  Matrix<ComplexScalar, 2, 2> t = m_matT.template block<2, 2>(iu - 1, iu - 1);
  const RealScalar normt = t.cwiseAbs().sum();
  const auto factors = internal::safe_scaling<RealScalar>::scale_to(t, t, normt);
 
  ComplexScalar b = t.coeff(0, 1) * t.coeff(1, 0);
  ComplexScalar c = t.coeff(0, 0) - t.coeff(1, 1);
  ComplexScalar disc = sqrt(c * c + RealScalar(4) * b);
  ComplexScalar det = t.coeff(0, 0) * t.coeff(1, 1) - b;
  ComplexScalar trace = t.coeff(0, 0) + t.coeff(1, 1);
  ComplexScalar eival1 = (trace + disc) / RealScalar(2);
  ComplexScalar eival2 = (trace - disc) / RealScalar(2);
  RealScalar eival1_norm = numext::norm1(eival1);
  RealScalar eival2_norm = numext::norm1(eival2);
  // A division by zero can only occur if eival1==eival2==0.
  // In this case, det==0, and all we have to do is checking that eival2_norm!=0
  if (eival1_norm > eival2_norm)
    eival2 = det / eival1;
  else if (!numext::is_exactly_zero(eival2_norm))
    eival1 = det / eival2;
 
  // choose the eigenvalue closest to the bottom entry of the diagonal
  ComplexScalar shift = numext::norm1(eival1 - t.coeff(1, 1)) < numext::norm1(eival2 - t.coeff(1, 1)) ? eival1 : eival2;
  internal::safe_scaling<RealScalar>::unscale_in_place(shift, factors);
  return shift;
}
 
template <typename MatrixType>
template <typename InputType>
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::compute(const EigenBase<InputType>& matrix, bool computeU) {
  m_matUisUptodate = false;
  eigen_assert(matrix.cols() == matrix.rows());
 
  if (matrix.cols() <= 1) {
    m_matT = matrix.derived().template cast<ComplexScalar>();
    if (computeU) m_matU = ComplexMatrixType::Identity(matrix.rows(), matrix.cols());
    m_info = Success;
    m_isInitialized = true;
    m_matUisUptodate = computeU;
    return *this;
  }
 
  // Reduce to Hessenberg form at unit scale, as RealSchur does: HessenbergDecomposition treats a subdiagonal tail
  // whose squared norm underflows as already zero, so an unscaled matrix near the bottom of the exponent range loses
  // its whole subdiagonal. For binary scalars the scale is the power of two just below the largest coefficient, clamped
  // so that it and its reciprocal are normal; when maxCoeff < min/eps, both steps round through integer significands,
  // which FTZ/DAZ and ARMv7 NEON cannot flush. A scale above one underflows coefficients more than the exponent range
  // below the largest one, a perturbation bounded by the smallest subnormal relative to that largest coefficient.
  // maxCoeff propagates NaN; a zero or non-finite maxCoeff carries no usable exponent and is left unscaled.
  const RealScalar maxCoeff = internal::safe_scaling<RealScalar>::recover_flushed_max_coeff(
      matrix.derived(), matrix.derived().cwiseAbs().template maxCoeff<PropagateNaN>());
  const internal::safe_scaling_factors<RealScalar> factors = internal::safe_scaling<RealScalar>::with_scaled(
      matrix.derived(), maxCoeff, [&](const auto& scaled) { m_hess.run(*this, scaled, computeU); });
  computeFromHessenberg(m_matT, m_matU, computeU);
  // m_matU is unitary either way; only the triangular factor carries the scale.
  internal::safe_scaling<RealScalar>::unscale_in_place(m_matT, maxCoeff, factors);
  return *this;
}

template <typename MatrixType>
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeInPlace(bool computeU) {
  m_matUisUptodate = false;
  eigen_assert(m_matT.cols() == m_matT.rows());
 
  if (m_matT.cols() <= 1) {
    if (computeU) m_matU = ComplexMatrixType::Identity(m_matT.rows(), m_matT.cols());
    m_info = Success;
    m_isInitialized = true;
    m_matUisUptodate = computeU;
    return *this;
  }
 
  const RealScalar maxCoeff = internal::safe_scaling<RealScalar>::recover_flushed_max_coeff(
      m_matT, m_matT.cwiseAbs().template maxCoeff<PropagateNaN>());
  const internal::safe_scaling_factors<RealScalar> factors =
      internal::safe_scaling<RealScalar>::scale_in_place(m_matT, maxCoeff);
  m_hess.runInPlace(*this, computeU);
  computeFromHessenberg(m_matT, m_matU, computeU);
  internal::safe_scaling<RealScalar>::unscale_in_place(m_matT, maxCoeff, factors);
  return *this;
}
 
template <typename MatrixType>
template <typename HessMatrixType, typename OrthMatrixType>
ComplexSchur<MatrixType>& ComplexSchur<MatrixType>::computeFromHessenberg(const HessMatrixType& matrixH,
                                                                          const OrthMatrixType& matrixQ,
                                                                          bool computeU) {
  if (!internal::is_same_dense(m_matT, matrixH)) m_matT = matrixH;
  if (computeU && !internal::is_same_dense(m_matU, matrixQ)) m_matU = matrixQ;
  reduceToTriangularForm(computeU);
  return *this;
}
namespace internal {
 
/* Reduce given matrix to Hessenberg form, and hold the workspace this takes. */
template <typename MatrixType, bool IsComplex>
struct complex_schur_reduce_to_hessenberg {
  // this is the implementation for the case IsComplex = true: the reduction runs in place in m_matT.
  using Schur = ComplexSchur<MatrixType>;
  using Scalar = typename Schur::ComplexScalar;
  using CoeffVectorType = Matrix<Scalar, Schur::RowsAtCompileTime == Dynamic ? Dynamic : Schur::RowsAtCompileTime - 1,
                                 1, Schur::Options & ~RowMajor,
                                 Schur::MaxRowsAtCompileTime == Dynamic ? Dynamic : Schur::MaxRowsAtCompileTime - 1, 1>;
  using WorkspaceType =
      Matrix<Scalar, Schur::ColsAtCompileTime, 1, Schur::Options & ~RowMajor, Schur::MaxColsAtCompileTime, 1>;
 
  explicit complex_schur_reduce_to_hessenberg(Index size) : m_workspace(size) {
    if (size > 1) m_hCoeffs.resize(size - 1);
  }
 
  template <typename InputType>
  void run(Schur& _this, const InputType& matrix, bool computeU) {
    _this.m_matT = matrix;
    runInPlace(_this, computeU);
  }
 
  void runInPlace(Schur& _this, bool computeU) {
    hessenberg_decomposition_inplace(_this.m_matT, m_hCoeffs, m_workspace, _this.m_matU, computeU);
  }
 
  CoeffVectorType m_hCoeffs;
  WorkspaceType m_workspace;
};
 
template <typename MatrixType>
struct complex_schur_reduce_to_hessenberg<MatrixType, false> {
  using Schur = ComplexSchur<MatrixType>;
 
  explicit complex_schur_reduce_to_hessenberg(Index size) : m_hess(size) {}
 
  template <typename InputType>
  void run(Schur& _this, const InputType& matrix, bool computeU) {
    using ComplexScalar = typename Schur::ComplexScalar;
 
    // Note: m_hess is over RealScalar; m_matT and m_matU is over ComplexScalar
    m_hess.compute(matrix);
    _this.m_matT = m_hess.matrixH().template cast<ComplexScalar>();
    if (computeU) {
      // TODO: this temporary allocation could potentially be avoided.
      MatrixType Q = m_hess.matrixQ();
      _this.m_matU = Q.template cast<ComplexScalar>();
    }
  }
 
  HessenbergDecomposition<MatrixType> m_hess;
};
 
}  // end namespace internal
 
// Reduce the Hessenberg matrix m_matT to triangular form by QR iteration.
template <typename MatrixType>
void ComplexSchur<MatrixType>::reduceToTriangularForm(bool computeU) {
  Index maxIters = m_maxIters;
  if (maxIters == -1) maxIters = m_maxIterationsPerRow * m_matT.rows();
 
  // The matrix m_matT is divided in three parts.
  // Rows 0,...,il-1 are decoupled from the rest because m_matT(il,il-1) is zero.
  // Rows il,...,iu is the part we are working on (the active submatrix).
  // Rows iu+1,...,end are already brought in triangular form.
  Index iu = m_matT.cols() - 1;
  Index il;
  Index iter = 0;       // number of iterations we are working on the (iu,iu) element
  Index totalIter = 0;  // number of iterations for whole matrix
 
  while (true) {
    // find iu, the bottom row of the active submatrix
    while (iu > 0) {
      if (!subdiagonalEntryIsNeglegible(iu - 1)) break;
      iter = 0;
      --iu;
    }
 
    // if iu is zero then we are done; the whole matrix is triangularized
    if (iu == 0) break;
 
    // if we spent too many iterations, we give up
    iter++;
    totalIter++;
    if (totalIter > maxIters) break;
 
    // find il, the top row of the active submatrix
    il = iu - 1;
    while (il > 0 && !subdiagonalEntryIsNeglegible(il - 1)) {
      --il;
    }
 
    /* perform the QR step using Givens rotations. The first rotation
       creates a bulge; the (il+2,il) element becomes nonzero. This
       bulge is chased down to the bottom of the active submatrix. */
 
    ComplexScalar shift = computeShift(iu, iter);
    JacobiRotation<ComplexScalar> rot;
    rot.makeGivens(m_matT.coeff(il, il) - shift, m_matT.coeff(il + 1, il));
    m_matT.rightCols(m_matT.cols() - il).applyOnTheLeft(il, il + 1, rot.adjoint());
    m_matT.topRows((std::min)(il + 2, iu) + 1).applyOnTheRight(il, il + 1, rot);
    if (computeU) m_matU.applyOnTheRight(il, il + 1, rot);
 
    for (Index i = il + 1; i < iu; i++) {
      rot.makeGivens(m_matT.coeffRef(i, i - 1), m_matT.coeffRef(i + 1, i - 1), &m_matT.coeffRef(i, i - 1));
      m_matT.coeffRef(i + 1, i - 1) = ComplexScalar(0);
      m_matT.rightCols(m_matT.cols() - i).applyOnTheLeft(i, i + 1, rot.adjoint());
      m_matT.topRows((std::min)(i + 2, iu) + 1).applyOnTheRight(i, i + 1, rot);
      if (computeU) m_matU.applyOnTheRight(i, i + 1, rot);
    }
  }
 
  if (totalIter <= maxIters)
    m_info = Success;
  else
    m_info = NoConvergence;
 
  m_isInitialized = true;
  m_matUisUptodate = computeU;
}
 
}  // end namespace Eigen
 
#endif  // EIGEN_COMPLEX_SCHUR_H