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 // Eigen is free software; you can redistribute it and/or
00008 // modify it under the terms of the GNU Lesser General Public
00009 // License as published by the Free Software Foundation; either
00010 // version 3 of the License, or (at your option) any later version.
00011 //
00012 // Alternatively, you can redistribute it and/or
00013 // modify it under the terms of the GNU General Public License as
00014 // published by the Free Software Foundation; either version 2 of
00015 // the License, or (at your option) any later version.
00016 //
00017 // Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
00018 // WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00019 // FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
00020 // GNU General Public License for more details.
00021 //
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License and a copy of the GNU General Public License along with
00024 // Eigen. If not, see <http://www.gnu.org/licenses/>.
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     // let's consider the vector v = i-th column starting at position i+1
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     // Apply similarity transformation to remaining columns,
00323     // i.e., compute A = H A H'
00324 
00325     // A = H A
00326     matA.bottomRightCorner(remainingSize, remainingSize)
00327         .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
00328 
00329     // A = A H'
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 } // end namespace internal
00385 
00386 } // end namespace Eigen
00387 
00388 #endif // EIGEN_HESSENBERGDECOMPOSITION_H