blob: fafab990ea46862e301a11a69438e5c872cfd07f [file] [log] [blame]
 // This file is part of Eigen, a lightweight C++ template library // for linear algebra. // // Copyright (C) 2008-2009 Gael Guennebaud // Copyright (C) 2010 Jitse Niesen // // 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_HESSENBERGDECOMPOSITION_H #define EIGEN_HESSENBERGDECOMPOSITION_H #include "./InternalHeaderCheck.h" namespace Eigen { namespace internal { template struct HessenbergDecompositionMatrixHReturnType; template struct traits > { typedef MatrixType ReturnType; }; } /** \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 class HessenbergDecomposition { public: /** \brief Synonym for the template parameter \p MatrixType_. */ typedef MatrixType_ MatrixType; enum { Size = MatrixType::RowsAtCompileTime, SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, Options = MatrixType::Options, MaxSize = MatrixType::MaxRowsAtCompileTime, MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 }; /** \brief Scalar type for matrices of type #MatrixType. */ typedef typename MatrixType::Scalar Scalar; typedef Eigen::Index 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-side * type. */ typedef Matrix CoeffVectorType; /** \brief Return type of matrixQ() */ typedef HouseholderSequence> HouseholderSequenceType; typedef internal::HessenbergDecompositionMatrixHReturnType MatrixHReturnType; /** \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 explicit HessenbergDecomposition(const EigenBase& matrix) : m_matrix(matrix.derived()), m_temp(matrix.rows()), m_isInitialized(false) { if(matrix.rows()<2) { m_isInitialized = true; return; } m_hCoeffs.resize(matrix.rows()-1,1); _compute(m_matrix, m_hCoeffs, m_temp); m_isInitialized = true; } /** \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, %Matrix * Computations). The cost is \f$10n^3/3 \f$ flops, where \f$n \f$ * denotes the size of the given matrix. * * This method reuses of the allocated data in the HessenbergDecomposition * object. * * Example: \include HessenbergDecomposition_compute.cpp * Output: \verbinclude HessenbergDecomposition_compute.out */ template HessenbergDecomposition& compute(const EigenBase& matrix) { m_matrix = matrix.derived(); if(matrix.rows()<2) { m_isInitialized = true; return *this; } m_hCoeffs.resize(matrix.rows()-1,1); _compute(m_matrix, m_hCoeffs, m_temp); m_isInitialized = true; 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: typedef Matrix VectorType; typedef typename NumTraits::Real RealScalar; static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); protected: MatrixType m_matrix; CoeffVectorType m_hCoeffs; VectorType m_temp; bool m_isInitialized; }; /** \internal * Performs a tridiagonal decomposition of \a matA in place. * * \param matA the input selfadjoint matrix * \param hCoeffs returned Householder coefficients * * The result is written in the lower triangular part of \a matA. * * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1. * * \sa packedMatrix() */ template void HessenbergDecomposition::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) { eigen_assert(matA.rows()==matA.cols()); Index n = matA.rows(); temp.resize(n); for (Index i = 0; i struct HessenbergDecompositionMatrixHReturnType : public ReturnByValue > { public: /** \brief Constructor. * * \param[in] hess Hessenberg decomposition */ HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition& 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 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().setZero(); } Index rows() const { return m_hess.packedMatrix().rows(); } Index cols() const { return m_hess.packedMatrix().cols(); } protected: const HessenbergDecomposition& m_hess; }; } // end namespace internal } // end namespace Eigen #endif // EIGEN_HESSENBERGDECOMPOSITION_H