00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 00005 // Copyright (C) 2010 Jitse Niesen <jitse@maths.leeds.ac.uk> 00006 // 00007 // Eigen is free software; you can redistribute it and/or 00008 // modify it under the terms of the GNU Lesser General Public 00009 // License as published by the Free Software Foundation; either 00010 // version 3 of the License, or (at your option) any later version. 00011 // 00012 // Alternatively, you can redistribute it and/or 00013 // modify it under the terms of the GNU General Public License as 00014 // published by the Free Software Foundation; either version 2 of 00015 // the License, or (at your option) any later version. 00016 // 00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY 00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the 00020 // GNU General Public License for more details. 00021 // 00022 // You should have received a copy of the GNU Lesser General Public 00023 // License and a copy of the GNU General Public License along with 00024 // Eigen. If not, see <http://www.gnu.org/licenses/>. 00025 00026 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H 00027 #define EIGEN_HESSENBERGDECOMPOSITION_H 00028 00029 namespace internal { 00030 00031 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType; 00032 template<typename MatrixType> 00033 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> > 00034 { 00035 typedef MatrixType ReturnType; 00036 }; 00037 00038 } 00039 00040 /** \eigenvalues_module \ingroup Eigenvalues_Module 00041 * 00042 * 00043 * \class HessenbergDecomposition 00044 * 00045 * \brief Reduces a square matrix to Hessenberg form by an orthogonal similarity transformation 00046 * 00047 * \tparam _MatrixType the type of the matrix of which we are computing the Hessenberg decomposition 00048 * 00049 * This class performs an Hessenberg decomposition of a matrix \f$ A \f$. In 00050 * the real case, the Hessenberg decomposition consists of an orthogonal 00051 * matrix \f$ Q \f$ and a Hessenberg matrix \f$ H \f$ such that \f$ A = Q H 00052 * Q^T \f$. An orthogonal matrix is a matrix whose inverse equals its 00053 * transpose (\f$ Q^{-1} = Q^T \f$). A Hessenberg matrix has zeros below the 00054 * subdiagonal, so it is almost upper triangular. The Hessenberg decomposition 00055 * of a complex matrix is \f$ A = Q H Q^* \f$ with \f$ Q \f$ unitary (that is, 00056 * \f$ Q^{-1} = Q^* \f$). 00057 * 00058 * Call the function compute() to compute the Hessenberg decomposition of a 00059 * given matrix. Alternatively, you can use the 00060 * HessenbergDecomposition(const MatrixType&) constructor which computes the 00061 * Hessenberg decomposition at construction time. Once the decomposition is 00062 * computed, you can use the matrixH() and matrixQ() functions to construct 00063 * the matrices H and Q in the decomposition. 00064 * 00065 * The documentation for matrixH() contains an example of the typical use of 00066 * this class. 00067 * 00068 * \sa class ComplexSchur, class Tridiagonalization, \ref QR_Module "QR Module" 00069 */ 00070 template<typename _MatrixType> class HessenbergDecomposition 00071 { 00072 public: 00073 00074 /** \brief Synonym for the template parameter \p _MatrixType. */ 00075 typedef _MatrixType MatrixType; 00076 00077 enum { 00078 Size = MatrixType::RowsAtCompileTime, 00079 SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, 00080 Options = MatrixType::Options, 00081 MaxSize = MatrixType::MaxRowsAtCompileTime, 00082 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 00083 }; 00084 00085 /** \brief Scalar type for matrices of type #MatrixType. */ 00086 typedef typename MatrixType::Scalar Scalar; 00087 typedef typename MatrixType::Index Index; 00088 00089 /** \brief Type for vector of Householder coefficients. 00090 * 00091 * This is column vector with entries of type #Scalar. The length of the 00092 * vector is one less than the size of #MatrixType, if it is a fixed-side 00093 * type. 00094 */ 00095 typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; 00096 00097 /** \brief Return type of matrixQ() */ 00098 typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType; 00099 00100 typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType; 00101 00102 /** \brief Default constructor; the decomposition will be computed later. 00103 * 00104 * \param [in] size The size of the matrix whose Hessenberg decomposition will be computed. 00105 * 00106 * The default constructor is useful in cases in which the user intends to 00107 * perform decompositions via compute(). The \p size parameter is only 00108 * used as a hint. It is not an error to give a wrong \p size, but it may 00109 * impair performance. 00110 * 00111 * \sa compute() for an example. 00112 */ 00113 HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) 00114 : m_matrix(size,size), 00115 m_temp(size), 00116 m_isInitialized(false) 00117 { 00118 if(size>1) 00119 m_hCoeffs.resize(size-1); 00120 } 00121 00122 /** \brief Constructor; computes Hessenberg decomposition of given matrix. 00123 * 00124 * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. 00125 * 00126 * This constructor calls compute() to compute the Hessenberg 00127 * decomposition. 00128 * 00129 * \sa matrixH() for an example. 00130 */ 00131 HessenbergDecomposition(const MatrixType& matrix) 00132 : m_matrix(matrix), 00133 m_temp(matrix.rows()), 00134 m_isInitialized(false) 00135 { 00136 if(matrix.rows()<2) 00137 { 00138 m_isInitialized = true; 00139 return; 00140 } 00141 m_hCoeffs.resize(matrix.rows()-1,1); 00142 _compute(m_matrix, m_hCoeffs, m_temp); 00143 m_isInitialized = true; 00144 } 00145 00146 /** \brief Computes Hessenberg decomposition of given matrix. 00147 * 00148 * \param[in] matrix Square matrix whose Hessenberg decomposition is to be computed. 00149 * \returns Reference to \c *this 00150 * 00151 * The Hessenberg decomposition is computed by bringing the columns of the 00152 * matrix successively in the required form using Householder reflections 00153 * (see, e.g., Algorithm 7.4.2 in Golub \& Van Loan, <i>%Matrix 00154 * Computations</i>). The cost is \f$ 10n^3/3 \f$ flops, where \f$ n \f$ 00155 * denotes the size of the given matrix. 00156 * 00157 * This method reuses of the allocated data in the HessenbergDecomposition 00158 * object. 00159 * 00160 * Example: \include HessenbergDecomposition_compute.cpp 00161 * Output: \verbinclude HessenbergDecomposition_compute.out 00162 */ 00163 HessenbergDecomposition& compute(const MatrixType& matrix) 00164 { 00165 m_matrix = matrix; 00166 if(matrix.rows()<2) 00167 { 00168 m_isInitialized = true; 00169 return *this; 00170 } 00171 m_hCoeffs.resize(matrix.rows()-1,1); 00172 _compute(m_matrix, m_hCoeffs, m_temp); 00173 m_isInitialized = true; 00174 return *this; 00175 } 00176 00177 /** \brief Returns the Householder coefficients. 00178 * 00179 * \returns a const reference to the vector of Householder coefficients 00180 * 00181 * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 00182 * or the member function compute(const MatrixType&) has been called 00183 * before to compute the Hessenberg decomposition of a matrix. 00184 * 00185 * The Householder coefficients allow the reconstruction of the matrix 00186 * \f$ Q \f$ in the Hessenberg decomposition from the packed data. 00187 * 00188 * \sa packedMatrix(), \ref Householder_Module "Householder module" 00189 */ 00190 const CoeffVectorType& householderCoefficients() const 00191 { 00192 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00193 return m_hCoeffs; 00194 } 00195 00196 /** \brief Returns the internal representation of the decomposition 00197 * 00198 * \returns a const reference to a matrix with the internal representation 00199 * of the decomposition. 00200 * 00201 * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 00202 * or the member function compute(const MatrixType&) has been called 00203 * before to compute the Hessenberg decomposition of a matrix. 00204 * 00205 * The returned matrix contains the following information: 00206 * - the upper part and lower sub-diagonal represent the Hessenberg matrix H 00207 * - the rest of the lower part contains the Householder vectors that, combined with 00208 * Householder coefficients returned by householderCoefficients(), 00209 * allows to reconstruct the matrix Q as 00210 * \f$ Q = H_{N-1} \ldots H_1 H_0 \f$. 00211 * Here, the matrices \f$ H_i \f$ are the Householder transformations 00212 * \f$ H_i = (I - h_i v_i v_i^T) \f$ 00213 * where \f$ h_i \f$ is the \f$ i \f$th Householder coefficient and 00214 * \f$ v_i \f$ is the Householder vector defined by 00215 * \f$ v_i = [ 0, \ldots, 0, 1, M(i+2,i), \ldots, M(N-1,i) ]^T \f$ 00216 * with M the matrix returned by this function. 00217 * 00218 * See LAPACK for further details on this packed storage. 00219 * 00220 * Example: \include HessenbergDecomposition_packedMatrix.cpp 00221 * Output: \verbinclude HessenbergDecomposition_packedMatrix.out 00222 * 00223 * \sa householderCoefficients() 00224 */ 00225 const MatrixType& packedMatrix() const 00226 { 00227 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00228 return m_matrix; 00229 } 00230 00231 /** \brief Reconstructs the orthogonal matrix Q in the decomposition 00232 * 00233 * \returns object representing the matrix Q 00234 * 00235 * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 00236 * or the member function compute(const MatrixType&) has been called 00237 * before to compute the Hessenberg decomposition of a matrix. 00238 * 00239 * This function returns a light-weight object of template class 00240 * HouseholderSequence. You can either apply it directly to a matrix or 00241 * you can convert it to a matrix of type #MatrixType. 00242 * 00243 * \sa matrixH() for an example, class HouseholderSequence 00244 */ 00245 HouseholderSequenceType matrixQ() const 00246 { 00247 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00248 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) 00249 .setLength(m_matrix.rows() - 1) 00250 .setShift(1); 00251 } 00252 00253 /** \brief Constructs the Hessenberg matrix H in the decomposition 00254 * 00255 * \returns expression object representing the matrix H 00256 * 00257 * \pre Either the constructor HessenbergDecomposition(const MatrixType&) 00258 * or the member function compute(const MatrixType&) has been called 00259 * before to compute the Hessenberg decomposition of a matrix. 00260 * 00261 * The object returned by this function constructs the Hessenberg matrix H 00262 * when it is assigned to a matrix or otherwise evaluated. The matrix H is 00263 * constructed from the packed matrix as returned by packedMatrix(): The 00264 * upper part (including the subdiagonal) of the packed matrix contains 00265 * the matrix H. It may sometimes be better to directly use the packed 00266 * matrix instead of constructing the matrix H. 00267 * 00268 * Example: \include HessenbergDecomposition_matrixH.cpp 00269 * Output: \verbinclude HessenbergDecomposition_matrixH.out 00270 * 00271 * \sa matrixQ(), packedMatrix() 00272 */ 00273 MatrixHReturnType matrixH() const 00274 { 00275 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00276 return MatrixHReturnType(*this); 00277 } 00278 00279 private: 00280 00281 typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType; 00282 typedef typename NumTraits<Scalar>::Real RealScalar; 00283 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); 00284 00285 protected: 00286 MatrixType m_matrix; 00287 CoeffVectorType m_hCoeffs; 00288 VectorType m_temp; 00289 bool m_isInitialized; 00290 }; 00291 00292 /** \internal 00293 * Performs a tridiagonal decomposition of \a matA in place. 00294 * 00295 * \param matA the input selfadjoint matrix 00296 * \param hCoeffs returned Householder coefficients 00297 * 00298 * The result is written in the lower triangular part of \a matA. 00299 * 00300 * Implemented from Golub's "%Matrix Computations", algorithm 8.3.1. 00301 * 00302 * \sa packedMatrix() 00303 */ 00304 template<typename MatrixType> 00305 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) 00306 { 00307 assert(matA.rows()==matA.cols()); 00308 Index n = matA.rows(); 00309 temp.resize(n); 00310 for (Index i = 0; i<n-1; ++i) 00311 { 00312 // let's consider the vector v = i-th column starting at position i+1 00313 Index remainingSize = n-i-1; 00314 RealScalar beta; 00315 Scalar h; 00316 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); 00317 matA.col(i).coeffRef(i+1) = beta; 00318 hCoeffs.coeffRef(i) = h; 00319 00320 // Apply similarity transformation to remaining columns, 00321 // i.e., compute A = H A H' 00322 00323 // A = H A 00324 matA.bottomRightCorner(remainingSize, remainingSize) 00325 .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0)); 00326 00327 // A = A H' 00328 matA.rightCols(remainingSize) 00329 .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0)); 00330 } 00331 } 00332 00333 namespace internal { 00334 00335 /** \eigenvalues_module \ingroup Eigenvalues_Module 00336 * 00337 * 00338 * \brief Expression type for return value of HessenbergDecomposition::matrixH() 00339 * 00340 * \tparam MatrixType type of matrix in the Hessenberg decomposition 00341 * 00342 * Objects of this type represent the Hessenberg matrix in the Hessenberg 00343 * decomposition of some matrix. The object holds a reference to the 00344 * HessenbergDecomposition class until the it is assigned or evaluated for 00345 * some other reason (the reference should remain valid during the life time 00346 * of this object). This class is the return type of 00347 * HessenbergDecomposition::matrixH(); there is probably no other use for this 00348 * class. 00349 */ 00350 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType 00351 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > 00352 { 00353 typedef typename MatrixType::Index Index; 00354 public: 00355 /** \brief Constructor. 00356 * 00357 * \param[in] hess Hessenberg decomposition 00358 */ 00359 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { } 00360 00361 /** \brief Hessenberg matrix in decomposition. 00362 * 00363 * \param[out] result Hessenberg matrix in decomposition \p hess which 00364 * was passed to the constructor 00365 */ 00366 template <typename ResultType> 00367 inline void evalTo(ResultType& result) const 00368 { 00369 result = m_hess.packedMatrix(); 00370 Index n = result.rows(); 00371 if (n>2) 00372 result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero(); 00373 } 00374 00375 Index rows() const { return m_hess.packedMatrix().rows(); } 00376 Index cols() const { return m_hess.packedMatrix().cols(); } 00377 00378 protected: 00379 const HessenbergDecomposition<MatrixType>& m_hess; 00380 }; 00381 00382 } 00383 00384 #endif // EIGEN_HESSENBERGDECOMPOSITION_H
Page generated by Doxygen 1.7.3 for MRPT 0.9.4 SVN: at Sat Mar 26 06:16:28 UTC 2011 |