// This file is part of Eigen, a lightweight C++ template library
// for linear algebra.
//
// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
// Copyright (C) 2010 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_HESSENBERGDECOMPOSITION_H
#define EIGEN_HESSENBERGDECOMPOSITION_H

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

namespace Eigen {

namespace internal {

template <typename MatrixType>
struct HessenbergDecompositionMatrixHReturnType;
template <typename MatrixType>
struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType>> {
  using ReturnType = typename MatrixType::PlainObject;
};

template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp);

}  // namespace internal

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 *
 * \class HessenbergDecomposition
 *
 * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation
 *
 * \tparam MatrixType_ the type of the matrix of which we are computing the Hessenberg decomposition
 *
 * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In
 * the real case, the Hessenberg decomposition consists of an orthogonal
 * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H
 * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its
 * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the
 * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition
 * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is,
 * \f$ Q^{-1} = Q^* \f$).
 *
 * Call the function compute() to compute the Hessenberg decomposition of a
 * given matrix. Alternatively, you can use the
 * HessenbergDecomposition(const MatrixType&) constructor which computes the
 * Hessenberg decomposition at construction time. Once the decomposition is
 * computed, you can use the matrixH() and matrixQ() functions to construct
 * the matrices H and Q in the decomposition.
 *
 * The documentation for matrixH() contains an example of the typical use of
 * this class.
 *
 * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module"
 */
template <typename MatrixType_>
class HessenbergDecomposition {
 public:
  /** \brief Synonym for the template parameter \p MatrixType_. */
  using MatrixType = MatrixType_;

  enum {
    Size = MatrixType::RowsAtCompileTime,
    SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
    Options = internal::plain_object_options<MatrixType>::value,
    MaxSize = MatrixType::MaxRowsAtCompileTime,
    MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
  };

  /** \brief Scalar type for matrices of type #MatrixType. */
  using Scalar = typename MatrixType::Scalar;
  using Index = Eigen::Index;  ///< \deprecated since Eigen 3.3

  /** \brief Type for vector of Householder coefficients.
   *
   * This is column vector with entries of type #Scalar. The length of the
   * vector is one less than the size of #MatrixType, if it is a fixed-size
   * type.
   */
  using CoeffVectorType = Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1>;

  /** \brief Return type of matrixQ() */
  using HouseholderSequenceType =
      HouseholderSequence<MatrixType, internal::remove_all_t<typename CoeffVectorType::ConjugateReturnType>>;

  using MatrixHReturnType = internal::HessenbergDecompositionMatrixHReturnType<MatrixType>;

  /** \brief Default constructor; the decomposition will be computed later.
   *
   * \param [in] size  The size of the matrix whose Hessenberg 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 HessenbergDecomposition(Index size = Size == Dynamic ? 2 : Size)
      : m_matrix(size, size), m_temp(size), m_isInitialized(false) {
    if (size > 1) m_hCoeffs.resize(size - 1);
  }

  /** \brief Constructor; computes Hessenberg decomposition of given matrix.
   *
   * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
   *
   * This constructor calls compute() to compute the Hessenberg
   * decomposition.
   *
   * \sa matrixH() for an example.
   */
  template <typename InputType>
  explicit HessenbergDecomposition(const EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) {
    computeInPlace();
  }

  /** \brief Constructor for \link InplaceDecomposition inplace decomposition \endlink
   *
   * \param[in,out]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
   *
   * When \p MatrixType is a Ref<>, the decomposition is computed within the memory of \p matrix, which then holds
   * the packed representation returned by packedMatrix(). Otherwise this constructor behaves like
   * HessenbergDecomposition(const EigenBase<InputType>&).
   */
  template <typename InputType>
  explicit HessenbergDecomposition(EigenBase<InputType>& matrix)
      : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) {
    computeInPlace();
  }

  /** \brief Computes Hessenberg decomposition of given matrix.
   *
   * \param[in]  matrix  Square matrix whose Hessenberg decomposition is to be computed.
   * \returns    Reference to \c *this
   *
   * The Hessenberg decomposition is computed by bringing the columns of the
   * matrix successively in the required form using Householder reflections
   * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix
   * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$
   * denotes the size of the given matrix.
   *
   * This method reuses the allocated data in the HessenbergDecomposition
   * object.
   *
   * Example: \include HessenbergDecomposition_compute.cpp
   * Output: \verbinclude HessenbergDecomposition_compute.out
   */
  template <typename InputType>
  HessenbergDecomposition& compute(const EigenBase<InputType>& matrix) {
    m_matrix = matrix.derived();
    computeInPlace();
    return *this;
  }

  /** \brief Returns the Householder coefficients.
   *
   * \returns a const reference to the vector of Householder coefficients
   *
   * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
   * or the member function compute(const MatrixType&) has been called
   * before to compute the Hessenberg decomposition of a matrix.
   *
   * The Householder coefficients allow the reconstruction of the matrix
   * \f$ Q \f$ in the Hessenberg decomposition from the packed data.
   *
   * \sa packedMatrix(), \ref Householder_Module "Householder module"
   */
  const CoeffVectorType& householderCoefficients() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return m_hCoeffs;
  }

  /** \brief Returns the internal representation of the decomposition
   *
   *	\returns a const reference to a matrix with the internal representation
   *	         of the decomposition.
   *
   * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
   * or the member function compute(const MatrixType&) has been called
   * before to compute the Hessenberg decomposition of a matrix.
   *
   * The returned matrix contains the following information:
   *  - the upper part and lower sub-diagonal represent the Hessenberg matrix H
   *  - the rest of the lower part contains the Householder vectors that, combined with
   *    Householder coefficients returned by householderCoefficients(),
   *    allows to reconstruct the matrix Q as
   *       \f$ Q = H_{N-1} \ldots H_1 H_0 \f$.
   *    Here, the matrices \f$ H_i \f$ are the Householder transformations
   *       \f$ H_i = (I - h_i v_i v_i^T) \f$
   *    where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and
   *    \f$ v_i \f$ is the Householder vector defined by
   *       \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$
   *    with M the matrix returned by this function.
   *
   * See LAPACK for further details on this packed storage.
   *
   * Example: \include HessenbergDecomposition_packedMatrix.cpp
   * Output: \verbinclude HessenbergDecomposition_packedMatrix.out
   *
   * \sa householderCoefficients()
   */
  const MatrixType& packedMatrix() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return m_matrix;
  }

  /** \brief Reconstructs the orthogonal matrix Q in the decomposition
   *
   * \returns object representing the matrix Q
   *
   * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
   * or the member function compute(const MatrixType&) has been called
   * before to compute the Hessenberg decomposition of a matrix.
   *
   * This function returns a light-weight object of template class
   * HouseholderSequence. You can either apply it directly to a matrix or
   * you can convert it to a matrix of type #MatrixType.
   *
   * \sa matrixH() for an example, class HouseholderSequence
   */
  HouseholderSequenceType matrixQ() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()).setLength(m_matrix.rows() - 1).setShift(1);
  }

  /** \brief Constructs the Hessenberg matrix H in the decomposition
   *
   * \returns expression object representing the matrix H
   *
   * \pre Either the constructor HessenbergDecomposition(const MatrixType&)
   * or the member function compute(const MatrixType&) has been called
   * before to compute the Hessenberg decomposition of a matrix.
   *
   * The object returned by this function constructs the Hessenberg matrix H
   * when it is assigned to a matrix or otherwise evaluated. The matrix H is
   * constructed from the packed matrix as returned by packedMatrix(): The
   * upper part (including the subdiagonal) of the packed matrix contains
   * the matrix H. It may sometimes be better to directly use the packed
   * matrix instead of constructing the matrix H.
   *
   * Example: \include HessenbergDecomposition_matrixH.cpp
   * Output: \verbinclude HessenbergDecomposition_matrixH.out
   *
   * \sa matrixQ(), packedMatrix()
   */
  MatrixHReturnType matrixH() const {
    eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
    return MatrixHReturnType(*this);
  }

 private:
  using VectorType = Matrix<Scalar, 1, Size, int(Options) | int(RowMajor), 1, MaxSize>;
  using RealScalar = typename NumTraits<Scalar>::Real;

  void computeInPlace() {
    if (m_matrix.rows() < 2) {
      m_isInitialized = true;
      return;
    }
    m_hCoeffs.resize(m_matrix.rows() - 1, 1);
    internal::hessenberg_decomposition_inplace(m_matrix, m_hCoeffs, m_temp);
    m_isInitialized = true;
  }

 protected:
  MatrixType m_matrix;
  CoeffVectorType m_hCoeffs;
  VectorType m_temp;
  bool m_isInitialized;
};

namespace internal {

/** \internal
 * Performs a Hessenberg decomposition of \a matA in place.
 *
 * \param matA the input square matrix
 * \param hCoeffs returned Householder coefficients
 * \param temp workspace of at least matA.rows() coefficients
 *
 * The result is written in the whole of \a matA: the upper part, including
 * the subdiagonal, holds the Hessenberg matrix H, while the part strictly
 * below the subdiagonal holds the Householder vectors.
 *
 * Implemented from Golub's "%Matrix Computations", algorithm 7.4.2.
 *
 * \sa HessenbergDecomposition::packedMatrix()
 */
template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp) {
  using Scalar = typename MatrixType::Scalar;
  using RealScalar = typename NumTraits<Scalar>::Real;
  EIGEN_INTERNAL_ASSERT_SQUARE(matA);
  Index n = matA.rows();
  temp.resize(n);
  for (Index i = 0; i < n - 1; ++i) {
    // let's consider the vector v = i-th column starting at position i+1
    Index remainingSize = n - i - 1;
    RealScalar beta;
    Scalar h;
    auto householder = matA.col(i).tail(remainingSize);
    const RealScalar tailSqNorm =
        remainingSize == 1 ? RealScalar(0) : householder.tail(remainingSize - 1).unwind().squaredNorm();
    const RealScalar tol = (std::numeric_limits<RealScalar>::min)();
    // Preserve negligible subdiagonal entries instead of rotating them into much larger matrix coefficients.  The
    // latter can erase small eigenvalues through cancellation even though the reflector itself is accurate.
    if (tailSqNorm <= tol && numext::abs2(numext::imag(householder.coeff(0))) <= tol) {
      h = Scalar(0);
      beta = numext::real(householder.coeff(0));
      householder.tail(remainingSize - 1).setZero();
    } else {
      householder.makeHouseholderInPlace(h, beta);
    }
    matA.col(i).coeffRef(i + 1) = beta;
    hCoeffs.coeffRef(i) = h;

    // Apply similarity transformation to remaining columns,
    // i.e., compute A = H A H'

    // A = H A
    matA.bottomRightCorner(remainingSize, remainingSize)
        .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize - 1), h, &temp.coeffRef(0));

    // A = A H'
    matA.rightCols(remainingSize)
        .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize - 1), numext::conj(h), &temp.coeffRef(0));
  }
}

/** \internal
 * Reduces \a matA in place to the Hessenberg matrix H of \f$ A = Q H Q^* \f$.
 *
 * \param matA the input square matrix, overwritten by H
 * \param hCoeffs returned Householder coefficients
 * \param temp workspace of at least matA.rows() coefficients
 * \param matQ returned unitary matrix Q, if \a computeQ is true
 * \param computeQ whether to compute Q
 *
 * Unlike the packed form above, \a matA holds only H afterwards: Q is formed from the Householder vectors below the
 * subdiagonal before they are cleared.
 */
template <typename MatrixType, typename CoeffVectorType, typename WorkspaceType, typename MatrixQType>
void hessenberg_decomposition_inplace(MatrixType& matA, CoeffVectorType& hCoeffs, WorkspaceType& temp,
                                      MatrixQType& matQ, bool computeQ) {
  using HouseholderSequenceType =
      HouseholderSequence<MatrixType, remove_all_t<typename CoeffVectorType::ConjugateReturnType>>;
  const Index n = matA.rows();
  hCoeffs.resize(n - 1);
  hessenberg_decomposition_inplace(matA, hCoeffs, temp);
  if (computeQ) {
    Index firstReflector = 0;
    while (firstReflector < n - 1 && numext::is_exactly_zero(hCoeffs.coeff(firstReflector))) ++firstReflector;
    if (firstReflector == n - 1) {
      // All tau_i = 0: Q = I, without allocating block Householder factors.
      matQ.setIdentity(n, n);
    } else {
      HouseholderSequenceType(matA, hCoeffs.conjugate()).setLength(n - 1).setShift(1).evalTo(matQ, temp);
    }
  }
  if (n > 2) matA.bottomLeftCorner(n - 2, n - 2).template triangularView<Lower>().setZero();
}

/** \eigenvalues_module \ingroup Eigenvalues_Module
 *
 *
 * \brief Expression type for return value of HessenbergDecomposition::matrixH()
 *
 * \tparam MatrixType type of matrix in the Hessenberg decomposition
 *
 * Objects of this type represent the Hessenberg matrix in the Hessenberg
 * decomposition of some matrix. The object holds a reference to the
 * HessenbergDecomposition class until it is assigned or evaluated for
 * some other reason (the reference should remain valid during the life time
 * of this object). This class is the return type of
 * HessenbergDecomposition::matrixH(); there is no other intended use for this
 * class.
 */
template <typename MatrixType>
struct HessenbergDecompositionMatrixHReturnType
    : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType>> {
 public:
  /** \brief Constructor.
   *
   * \param[in] hess  Hessenberg decomposition
   */
  HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) {}

  /** \brief Hessenberg matrix in decomposition.
   *
   * \param[out] result  Hessenberg matrix in decomposition \p hess which
   *                     was passed to the constructor
   */
  template <typename ResultType>
  inline void evalTo(ResultType& result) const {
    result = m_hess.packedMatrix();
    Index n = result.rows();
    if (n > 2) result.bottomLeftCorner(n - 2, n - 2).template triangularView<Lower>().setZero();
  }

  Index rows() const { return m_hess.packedMatrix().rows(); }
  Index cols() const { return m_hess.packedMatrix().cols(); }

 protected:
  const HessenbergDecomposition<MatrixType>& m_hess;
};

}  // end namespace internal

}  // end namespace Eigen

#endif  // EIGEN_HESSENBERGDECOMPOSITION_H