$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 typename MatrixType::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 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 HessenbergDecomposition(const MatrixType& matrix) 00119 : m_matrix(matrix), 00120 m_temp(matrix.rows()), 00121 m_isInitialized(false) 00122 { 00123 if(matrix.rows()<2) 00124 { 00125 m_isInitialized = true; 00126 return; 00127 } 00128 m_hCoeffs.resize(matrix.rows()-1,1); 00129 _compute(m_matrix, m_hCoeffs, m_temp); 00130 m_isInitialized = true; 00131 } 00132 00150 HessenbergDecomposition& compute(const MatrixType& matrix) 00151 { 00152 m_matrix = matrix; 00153 if(matrix.rows()<2) 00154 { 00155 m_isInitialized = true; 00156 return *this; 00157 } 00158 m_hCoeffs.resize(matrix.rows()-1,1); 00159 _compute(m_matrix, m_hCoeffs, m_temp); 00160 m_isInitialized = true; 00161 return *this; 00162 } 00163 00177 const CoeffVectorType& householderCoefficients() const 00178 { 00179 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00180 return m_hCoeffs; 00181 } 00182 00212 const MatrixType& packedMatrix() const 00213 { 00214 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00215 return m_matrix; 00216 } 00217 00232 HouseholderSequenceType matrixQ() const 00233 { 00234 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00235 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate()) 00236 .setLength(m_matrix.rows() - 1) 00237 .setShift(1); 00238 } 00239 00260 MatrixHReturnType matrixH() const 00261 { 00262 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized."); 00263 return MatrixHReturnType(*this); 00264 } 00265 00266 private: 00267 00268 typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType; 00269 typedef typename NumTraits<Scalar>::Real RealScalar; 00270 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp); 00271 00272 protected: 00273 MatrixType m_matrix; 00274 CoeffVectorType m_hCoeffs; 00275 VectorType m_temp; 00276 bool m_isInitialized; 00277 }; 00278 00291 template<typename MatrixType> 00292 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp) 00293 { 00294 eigen_assert(matA.rows()==matA.cols()); 00295 Index n = matA.rows(); 00296 temp.resize(n); 00297 for (Index i = 0; i<n-1; ++i) 00298 { 00299 // let's consider the vector v = i-th column starting at position i+1 00300 Index remainingSize = n-i-1; 00301 RealScalar beta; 00302 Scalar h; 00303 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta); 00304 matA.col(i).coeffRef(i+1) = beta; 00305 hCoeffs.coeffRef(i) = h; 00306 00307 // Apply similarity transformation to remaining columns, 00308 // i.e., compute A = H A H' 00309 00310 // A = H A 00311 matA.bottomRightCorner(remainingSize, remainingSize) 00312 .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0)); 00313 00314 // A = A H' 00315 matA.rightCols(remainingSize) 00316 .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), numext::conj(h), &temp.coeffRef(0)); 00317 } 00318 } 00319 00320 namespace internal { 00321 00337 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType 00338 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> > 00339 { 00340 typedef typename MatrixType::Index Index; 00341 public: 00346 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { } 00347 00353 template <typename ResultType> 00354 inline void evalTo(ResultType& result) const 00355 { 00356 result = m_hess.packedMatrix(); 00357 Index n = result.rows(); 00358 if (n>2) 00359 result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero(); 00360 } 00361 00362 Index rows() const { return m_hess.packedMatrix().rows(); } 00363 Index cols() const { return m_hess.packedMatrix().cols(); } 00364 00365 protected: 00366 const HessenbergDecomposition<MatrixType>& m_hess; 00367 }; 00368 00369 } // end namespace internal 00370 00371 } // end namespace Eigen 00372 00373 #endif // EIGEN_HESSENBERGDECOMPOSITION_H