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