HessenbergDecomposition.h
Go to the documentation of this file.
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


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 11:58:27