![]() |
Eigen
3.3.3
|
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 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H 00012 #define EIGEN_HESSENBERGDECOMPOSITION_H 00013 00014 namespace Eigen { 00015 00016 namespace internal { 00017 00018 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType; 00019 template<typename MatrixType> 00020 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> > 00021 { 00022 typedef MatrixType ReturnType; 00023 }; 00024 00025 } 00026 00057 template<typename _MatrixType> class HessenbergDecomposition 00058 { 00059 public: 00060 00062 typedef _MatrixType MatrixType; 00063 00064 enum { 00065 Size = MatrixType::RowsAtCompileTime, 00066 SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1, 00067 Options = MatrixType::Options, 00068 MaxSize = MatrixType::MaxRowsAtCompileTime, 00069 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1 00070 }; 00071 00073 typedef typename MatrixType::Scalar Scalar; 00074 typedef Eigen::Index Index; 00075 00082 typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType; 00083 00085 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename CoeffVectorType::ConjugateReturnType>::type> HouseholderSequenceType; 00086 00087 typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType; 00088 00100 explicit HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size) 00101 : m_matrix(size,size), 00102 m_temp(size), 00103 m_isInitialized(false) 00104 { 00105 if(size>1) 00106 m_hCoeffs.resize(size-1); 00107 } 00108 00118 template<typename InputType> 00119 explicit HessenbergDecomposition(const EigenBase<InputType>& matrix) 00120 : m_matrix(matrix.derived()), 00121 m_temp(matrix.rows()), 00122 m_isInitialized(false) 00123 { 00124 if(matrix.rows()<2) 00125 { 00126 m_isInitialized = true; 00127 return; 00128 } 00129 m_hCoeffs.resize(matrix.rows()-1,1); 00130 _compute(m_matrix, m_hCoeffs, m_temp); 00131 m_isInitialized = true; 00132 } 00133 00151 template<typename InputType> 00152 HessenbergDecomposition& compute(const EigenBase<InputType>& matrix) 00153 { 00154 m_matrix = matrix.derived(); 00155 if(matrix.rows()<2) 00156 { 00157 m_isInitialized = true; 00158 return *this; 00159 } 00160 m_hCoeffs.resize(matrix.rows()-1,1); 00161 _compute(m_matrix, m_hCoeffs, m_temp); 00162 m_isInitialized = true; 00163 return *this; 00164 } 00165 00179 const CoeffVectorType& householderCoefficients() const 00180 { 00181 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00182 return m_hCoeffs; 00183 } 00184 00214 const MatrixType& packedMatrix() const 00215 { 00216 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00217 return m_matrix; 00218 } 00219 00234 HouseholderSequenceType matrixQ() const 00235 { 00236 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00237 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) 00238 .setLength(m_matrix.rows() - 1) 00239 .setShift(1); 00240 } 00241 00262 MatrixHReturnType matrixH() const 00263 { 00264 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00265 return MatrixHReturnType(*this); 00266 } 00267 00268 private: 00269 00270 typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType; 00271 typedef typename NumTraits<Scalar>::Real RealScalar; 00272 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); 00273 00274 protected: 00275 MatrixType m_matrix; 00276 CoeffVectorType m_hCoeffs; 00277 VectorType m_temp; 00278 bool m_isInitialized; 00279 }; 00280 00293 template<typename MatrixType> 00294 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) 00295 { 00296 eigen_assert(matA.rows()==matA.cols()); 00297 Index n = matA.rows(); 00298 temp.resize(n); 00299 for (Index i = 0; i<n-1; ++i) 00300 { 00301 // let's consider the vector v = i-th column starting at position i+1 00302 Index remainingSize = n-i-1; 00303 RealScalar beta; 00304 Scalar h; 00305 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); 00306 matA.col(i).coeffRef(i+1) = beta; 00307 hCoeffs.coeffRef(i) = h; 00308 00309 // Apply similarity transformation to remaining columns, 00310 // i.e., compute A = H A H' 00311 00312 // A = H A 00313 matA.bottomRightCorner(remainingSize, remainingSize) 00314 .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0)); 00315 00316 // A = A H' 00317 matA.rightCols(remainingSize) 00318 .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0)); 00319 } 00320 } 00321 00322 namespace internal { 00323 00339 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType 00340 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > 00341 { 00342 public: 00347 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { } 00348 00354 template <typename ResultType> 00355 inline void evalTo(ResultType& result) const 00356 { 00357 result = m_hess.packedMatrix(); 00358 Index n = result.rows(); 00359 if (n>2) 00360 result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero(); 00361 } 00362 00363 Index rows() const { return m_hess.packedMatrix().rows(); } 00364 Index cols() const { return m_hess.packedMatrix().cols(); } 00365 00366 protected: 00367 const HessenbergDecomposition<MatrixType>& m_hess; 00368 }; 00369 00370 } // end namespace internal 00371 00372 } // end namespace Eigen 00373 00374 #endif // EIGEN_HESSENBERGDECOMPOSITION_H