1c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This file is part of Eigen, a lightweight C++ template library 2c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// for linear algebra. 3c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 4c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 6c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// 7c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// This Source Code Form is subject to the terms of the Mozilla 8c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// Public License v. 2.0. If a copy of the MPL was not distributed 9c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath// with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 11c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#ifndef EIGEN_HESSENBERGDECOMPOSITION_H 12c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#define EIGEN_HESSENBERGDECOMPOSITION_H 13c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 14c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace Eigen { 15c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 16c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 17c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 18c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType; 19c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 20c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathstruct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> > 21c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 22c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef MatrixType ReturnType; 23c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 24c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 25c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 26c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 27c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 28c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 29c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 30c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \class HessenbergDecomposition 31c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 32c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation 33c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 34c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition 35c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 36c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In 37c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the real case, the Hessenberg decomposition consists of an orthogonal 38c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H 39c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its 40c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the 41c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition 42c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is, 43c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q^{-1} = Q^* \f$). 44c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 45c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Call the function compute() to compute the Hessenberg decomposition of a 46c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * given matrix. Alternatively, you can use the 47c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HessenbergDecomposition(const MatrixType&) constructor which computes the 48c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Hessenberg decomposition at construction time. Once the decomposition is 49c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * computed, you can use the matrixH() and matrixQ() functions to construct 50c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the matrices H and Q in the decomposition. 51c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 52c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The documentation for matrixH() contains an example of the typical use of 53c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * this class. 54c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 55c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module" 56c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 57c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename _MatrixType> class HessenbergDecomposition 58c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 59c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 60c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 61c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Synonym for the template parameter \p _MatrixType. */ 62c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef _MatrixType MatrixType; 63c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 64c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath enum { 65c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Size = MatrixType::RowsAtCompileTime, 66c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, 67c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Options = MatrixType::Options, 68c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxSize = MatrixType::MaxRowsAtCompileTime, 69c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 70c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath }; 71c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 72c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Scalar type for matrices of type #MatrixType. */ 73c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Scalar Scalar; 74c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 75c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 76c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Type for vector of Householder coefficients. 77c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 78c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This is column vector with entries of type #Scalar. The length of the 79c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * vector is one less than the size of #MatrixType, if it is a fixed-side 80c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * type. 81c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 82c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; 83c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 84c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Return type of matrixQ() */ 857faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType; 86c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 87c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType; 88c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 89c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Default constructor; the decomposition will be computed later. 90c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 91c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed. 92c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 93c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The default constructor is useful in cases in which the user intends to 94c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * perform decompositions via compute(). The \p size parameter is only 95c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * used as a hint. It is not an error to give a wrong \p size, but it may 96c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * impair performance. 97c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 98c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa compute() for an example. 99c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 100c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) 101c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(size,size), 102c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temp(size), 103c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 104c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 105c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(size>1) 106c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs.resize(size-1); 107c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 108c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 109c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor; computes Hessenberg decomposition of given matrix. 110c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 111c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. 112c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 113c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This constructor calls compute() to compute the Hessenberg 114c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition. 115c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 116c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa matrixH() for an example. 117c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 118c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecomposition(const MatrixType& matrix) 119c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath : m_matrix(matrix), 120c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_temp(matrix.rows()), 121c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized(false) 122c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 123c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.rows()<2) 124c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 125c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 126c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return; 127c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 128c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs.resize(matrix.rows()-1,1); 129c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _compute(m_matrix, m_hCoeffs, m_temp); 130c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 131c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 132c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 133c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Computes Hessenberg decomposition of given matrix. 134c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 135c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. 136c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns Reference to \c *this 137c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 138c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Hessenberg decomposition is computed by bringing the columns of the 139c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix successively in the required form using Householder reflections 140c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix 141c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$ 142c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * denotes the size of the given matrix. 143c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 144c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This method reuses of the allocated data in the HessenbergDecomposition 145c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * object. 146c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 147c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include HessenbergDecomposition_compute.cpp 148c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude HessenbergDecomposition_compute.out 149c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 150c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecomposition& compute(const MatrixType& matrix) 151c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 152c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_matrix = matrix; 153c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if(matrix.rows()<2) 154c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 155c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 156c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 157c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 158c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_hCoeffs.resize(matrix.rows()-1,1); 159c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath _compute(m_matrix, m_hCoeffs, m_temp); 160c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath m_isInitialized = true; 161c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return *this; 162c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 163c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 164c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the Householder coefficients. 165c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 166c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a const reference to the vector of Householder coefficients 167c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 168c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 169c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or the member function compute(const MatrixType&) has been called 170c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * before to compute the Hessenberg decomposition of a matrix. 171c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 172c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The Householder coefficients allow the reconstruction of the matrix 173c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q \f$ in the Hessenberg decomposition from the packed data. 174c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 175c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa packedMatrix(), \ref Householder_Module "Householder module" 176c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 177c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const CoeffVectorType& householderCoefficients() const 178c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 179c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 180c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_hCoeffs; 181c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 182c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 183c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Returns the internal representation of the decomposition 184c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 185c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns a const reference to a matrix with the internal representation 186c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of the decomposition. 187c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 188c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 189c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or the member function compute(const MatrixType&) has been called 190c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * before to compute the Hessenberg decomposition of a matrix. 191c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 192c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The returned matrix contains the following information: 193c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - the upper part and lower sub-diagonal represent the Hessenberg matrix H 194c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * - the rest of the lower part contains the Householder vectors that, combined with 195c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Householder coefficients returned by householderCoefficients(), 196c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * allows to reconstruct the matrix Q as 197c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. 198c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Here, the matrices \f$ H_i \f$ are the Householder transformations 199c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ H_i = (I - h_i v_i v_i^T) \f$ 200c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and 201c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i \f$ is the Householder vector defined by 202c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$ 203c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * with M the matrix returned by this function. 204c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 205c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * See LAPACK for further details on this packed storage. 206c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 207c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include HessenbergDecomposition_packedMatrix.cpp 208c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude HessenbergDecomposition_packedMatrix.out 209c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 210c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa householderCoefficients() 211c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 212c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const MatrixType& packedMatrix() const 213c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 214c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 215c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return m_matrix; 216c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 217c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 218c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Reconstructs the orthogonal matrix Q in the decomposition 219c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 220c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns object representing the matrix Q 221c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 222c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 223c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or the member function compute(const MatrixType&) has been called 224c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * before to compute the Hessenberg decomposition of a matrix. 225c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 226c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * This function returns a light-weight object of template class 227c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HouseholderSequence. You can either apply it directly to a matrix or 228c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * you can convert it to a matrix of type #MatrixType. 229c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 230c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa matrixH() for an example, class HouseholderSequence 231c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 232c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HouseholderSequenceType matrixQ() const 233c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 234c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 235c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) 236c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setLength(m_matrix.rows() - 1) 237c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .setShift(1); 238c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 239c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 240c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructs the Hessenberg matrix H in the decomposition 241c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 242c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \returns expression object representing the matrix H 243c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 244c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 245c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * or the member function compute(const MatrixType&) has been called 246c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * before to compute the Hessenberg decomposition of a matrix. 247c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 248c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The object returned by this function constructs the Hessenberg matrix H 249c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * when it is assigned to a matrix or otherwise evaluated. The matrix H is 250c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * constructed from the packed matrix as returned by packedMatrix(): The 251c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * upper part (including the subdiagonal) of the packed matrix contains 252c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * the matrix H. It may sometimes be better to directly use the packed 253c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * matrix instead of constructing the matrix H. 254c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 255c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Example: \include HessenbergDecomposition_matrixH.cpp 256c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Output: \verbinclude HessenbergDecomposition_matrixH.out 257c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 258c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa matrixQ(), packedMatrix() 259c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 260c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixHReturnType matrixH() const 261c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 262c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 263c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath return MatrixHReturnType(*this); 264c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 265c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 266c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath private: 267c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 268c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType; 269c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename NumTraits<Scalar>::Real RealScalar; 270c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); 271c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 272c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 273c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath MatrixType m_matrix; 274c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath CoeffVectorType m_hCoeffs; 275c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath VectorType m_temp; 276c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath bool m_isInitialized; 277c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 278c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 279c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \internal 280c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Performs a tridiagonal decomposition of \a matA in place. 281c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 282c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param matA the input selfadjoint matrix 283c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param hCoeffs returned Householder coefficients 284c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 285c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * The result is written in the lower triangular part of \a matA. 286c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 287c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1. 288c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 289c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \sa packedMatrix() 290c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 291c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> 292c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathvoid HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) 293c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 2947faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez eigen_assert(matA.rows()==matA.cols()); 295c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = matA.rows(); 296c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath temp.resize(n); 297c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath for (Index i = 0; i<n-1; ++i) 298c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 299c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // let's consider the vector v = i-th column starting at position i+1 300c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index remainingSize = n-i-1; 301c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath RealScalar beta; 302c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Scalar h; 303c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); 304c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.col(i).coeffRef(i+1) = beta; 305c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath hCoeffs.coeffRef(i) = h; 306c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 307c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // Apply similarity transformation to remaining columns, 308c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // i.e., compute A = H A H' 309c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 310c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A = H A 311c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.bottomRightCorner(remainingSize, remainingSize) 312c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0)); 313c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 314c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath // A = A H' 315c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath matA.rightCols(remainingSize) 3167faaa9f3f0df9d23790277834d426c3d992ac3baCarlos Hernandez .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0)); 317c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 318c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} 319c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 320c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathnamespace internal { 321c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 322c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath/** \eigenvalues_module \ingroup Eigenvalues_Module 323c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 324c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 325c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \brief Expression type for return value of HessenbergDecomposition::matrixH() 326c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 327c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \tparam MatrixType type of matrix in the Hessenberg decomposition 328c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 329c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * Objects of this type represent the Hessenberg matrix in the Hessenberg 330c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * decomposition of some matrix. The object holds a reference to the 331c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HessenbergDecomposition class until the it is assigned or evaluated for 332c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * some other reason (the reference should remain valid during the life time 333c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * of this object). This class is the return type of 334c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * HessenbergDecomposition::matrixH(); there is probably no other use for this 335c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * class. 336c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 337c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamathtemplate<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType 338c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath: public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > 339c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath{ 340c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath typedef typename MatrixType::Index Index; 341c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath public: 342c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Constructor. 343c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 344c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[in] hess Hessenberg decomposition 345c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 346c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { } 347c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 348c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath /** \brief Hessenberg matrix in decomposition. 349c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * 350c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * \param[out] result Hessenberg matrix in decomposition \p hess which 351c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath * was passed to the constructor 352c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath */ 353c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath template <typename ResultType> 354c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath inline void evalTo(ResultType& result) const 355c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath { 356c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result = m_hess.packedMatrix(); 357c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index n = result.rows(); 358c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath if (n>2) 359c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero(); 360c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath } 361c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 362c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index rows() const { return m_hess.packedMatrix().rows(); } 363c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath Index cols() const { return m_hess.packedMatrix().cols(); } 364c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 365c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath protected: 366c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath const HessenbergDecomposition<MatrixType>& m_hess; 367c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath}; 368c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 369c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace internal 370c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 371c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath} // end namespace Eigen 372c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath 373c981c48f5bc9aefeffc0bcb0cc3934c2fae179ddNarayan Kamath#endif // EIGEN_HESSENBERGDECOMPOSITION_H 374