Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026 #ifndef EIGEN_HESSENBERGDECOMPOSITION_H
00027 #define EIGEN_HESSENBERGDECOMPOSITION_H
00028
00029 namespace Eigen {
00030
00031 namespace internal {
00032
00033 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
00034 template<typename MatrixType>
00035 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00036 {
00037 typedef MatrixType ReturnType;
00038 };
00039
00040 }
00041
00072 template<typename _MatrixType> class HessenbergDecomposition
00073 {
00074 public:
00075
00077 typedef _MatrixType MatrixType;
00078
00079 enum {
00080 Size = MatrixType::RowsAtCompileTime,
00081 SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
00082 Options = MatrixType::Options,
00083 MaxSize = MatrixType::MaxRowsAtCompileTime,
00084 MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
00085 };
00086
00088 typedef typename MatrixType::Scalar Scalar;
00089 typedef typename MatrixType::Index Index;
00090
00097 typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00098
00100 typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00101
00102 typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
00103
00115 HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
00116 : m_matrix(size,size),
00117 m_temp(size),
00118 m_isInitialized(false)
00119 {
00120 if(size>1)
00121 m_hCoeffs.resize(size-1);
00122 }
00123
00133 HessenbergDecomposition(const MatrixType& matrix)
00134 : m_matrix(matrix),
00135 m_temp(matrix.rows()),
00136 m_isInitialized(false)
00137 {
00138 if(matrix.rows()<2)
00139 {
00140 m_isInitialized = true;
00141 return;
00142 }
00143 m_hCoeffs.resize(matrix.rows()-1,1);
00144 _compute(m_matrix, m_hCoeffs, m_temp);
00145 m_isInitialized = true;
00146 }
00147
00165 HessenbergDecomposition& compute(const MatrixType& matrix)
00166 {
00167 m_matrix = matrix;
00168 if(matrix.rows()<2)
00169 {
00170 m_isInitialized = true;
00171 return *this;
00172 }
00173 m_hCoeffs.resize(matrix.rows()-1,1);
00174 _compute(m_matrix, m_hCoeffs, m_temp);
00175 m_isInitialized = true;
00176 return *this;
00177 }
00178
00192 const CoeffVectorType& householderCoefficients() const
00193 {
00194 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00195 return m_hCoeffs;
00196 }
00197
00227 const MatrixType& packedMatrix() const
00228 {
00229 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00230 return m_matrix;
00231 }
00232
00247 HouseholderSequenceType matrixQ() const
00248 {
00249 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00250 return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00251 .setLength(m_matrix.rows() - 1)
00252 .setShift(1);
00253 }
00254
00275 MatrixHReturnType matrixH() const
00276 {
00277 eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00278 return MatrixHReturnType(*this);
00279 }
00280
00281 private:
00282
00283 typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
00284 typedef typename NumTraits<Scalar>::Real RealScalar;
00285 static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
00286
00287 protected:
00288 MatrixType m_matrix;
00289 CoeffVectorType m_hCoeffs;
00290 VectorType m_temp;
00291 bool m_isInitialized;
00292 };
00293
00306 template<typename MatrixType>
00307 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
00308 {
00309 assert(matA.rows()==matA.cols());
00310 Index n = matA.rows();
00311 temp.resize(n);
00312 for (Index i = 0; i<n-1; ++i)
00313 {
00314
00315 Index remainingSize = n-i-1;
00316 RealScalar beta;
00317 Scalar h;
00318 matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00319 matA.col(i).coeffRef(i+1) = beta;
00320 hCoeffs.coeffRef(i) = h;
00321
00322
00323
00324
00325
00326 matA.bottomRightCorner(remainingSize, remainingSize)
00327 .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
00328
00329
00330 matA.rightCols(remainingSize)
00331 .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
00332 }
00333 }
00334
00335 namespace internal {
00336
00352 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
00353 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00354 {
00355 typedef typename MatrixType::Index Index;
00356 public:
00361 HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
00362
00368 template <typename ResultType>
00369 inline void evalTo(ResultType& result) const
00370 {
00371 result = m_hess.packedMatrix();
00372 Index n = result.rows();
00373 if (n>2)
00374 result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
00375 }
00376
00377 Index rows() const { return m_hess.packedMatrix().rows(); }
00378 Index cols() const { return m_hess.packedMatrix().cols(); }
00379
00380 protected:
00381 const HessenbergDecomposition<MatrixType>& m_hess;
00382 };
00383
00384 }
00385
00386 }
00387
00388 #endif // EIGEN_HESSENBERGDECOMPOSITION_H