Eigen  3.3.3
HessenbergDecomposition.h
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
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends