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 internal {
00030   
00031 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType;
00032 template<typename MatrixType>
00033 struct traits<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00034 {
00035   typedef MatrixType ReturnType;
00036 };
00037 
00038 }
00039 
00070 template<typename _MatrixType> class HessenbergDecomposition
00071 {
00072   public:
00073 
00075     typedef _MatrixType MatrixType;
00076 
00077     enum {
00078       Size = MatrixType::RowsAtCompileTime,
00079       SizeMinusOne = Size == Dynamic ? Dynamic : Size - 1,
00080       Options = MatrixType::Options,
00081       MaxSize = MatrixType::MaxRowsAtCompileTime,
00082       MaxSizeMinusOne = MaxSize == Dynamic ? Dynamic : MaxSize - 1
00083     };
00084 
00086     typedef typename MatrixType::Scalar Scalar;
00087     typedef typename MatrixType::Index Index;
00088 
00095     typedef Matrix<Scalar, SizeMinusOne, 1, Options & ~RowMajor, MaxSizeMinusOne, 1> CoeffVectorType;
00096 
00098     typedef typename HouseholderSequence<MatrixType,CoeffVectorType>::ConjugateReturnType HouseholderSequenceType;
00099     
00100     typedef internal::HessenbergDecompositionMatrixHReturnType<MatrixType> MatrixHReturnType;
00101 
00113     HessenbergDecomposition(Index size = Size==Dynamic ? 2 : Size)
00114       : m_matrix(size,size),
00115         m_temp(size),
00116         m_isInitialized(false)
00117     {
00118       if(size>1)
00119         m_hCoeffs.resize(size-1);
00120     }
00121 
00131     HessenbergDecomposition(const MatrixType& matrix)
00132       : m_matrix(matrix),
00133         m_temp(matrix.rows()),
00134         m_isInitialized(false)
00135     {
00136       if(matrix.rows()<2)
00137       {
00138         m_isInitialized = true;
00139         return;
00140       }
00141       m_hCoeffs.resize(matrix.rows()-1,1);
00142       _compute(m_matrix, m_hCoeffs, m_temp);
00143       m_isInitialized = true;
00144     }
00145 
00163     HessenbergDecomposition& compute(const MatrixType& matrix)
00164     {
00165       m_matrix = matrix;
00166       if(matrix.rows()<2)
00167       {
00168         m_isInitialized = true;
00169         return *this;
00170       }
00171       m_hCoeffs.resize(matrix.rows()-1,1);
00172       _compute(m_matrix, m_hCoeffs, m_temp);
00173       m_isInitialized = true;
00174       return *this;
00175     }
00176 
00190     const CoeffVectorType& householderCoefficients() const
00191     {
00192       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00193       return m_hCoeffs;
00194     }
00195 
00225     const MatrixType& packedMatrix() const
00226     {
00227       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00228       return m_matrix;
00229     }
00230 
00245     HouseholderSequenceType matrixQ() const
00246     {
00247       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00248       return HouseholderSequenceType(m_matrix, m_hCoeffs.conjugate())
00249              .setLength(m_matrix.rows() - 1)
00250              .setShift(1);
00251     }
00252 
00273     MatrixHReturnType matrixH() const
00274     {
00275       eigen_assert(m_isInitialized && "HessenbergDecomposition is not initialized.");
00276       return MatrixHReturnType(*this);
00277     }
00278 
00279   private:
00280 
00281     typedef Matrix<Scalar, 1, Size, Options | RowMajor, 1, MaxSize> VectorType;
00282     typedef typename NumTraits<Scalar>::Real RealScalar;
00283     static void _compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp);
00284 
00285   protected:
00286     MatrixType m_matrix;
00287     CoeffVectorType m_hCoeffs;
00288     VectorType m_temp;
00289     bool m_isInitialized;
00290 };
00291 
00304 template<typename MatrixType>
00305 void HessenbergDecomposition<MatrixType>::_compute(MatrixType& matA, CoeffVectorType& hCoeffs, VectorType& temp)
00306 {
00307   assert(matA.rows()==matA.cols());
00308   Index n = matA.rows();
00309   temp.resize(n);
00310   for (Index i = 0; i<n-1; ++i)
00311   {
00312     // let's consider the vector v = i-th column starting at position i+1
00313     Index remainingSize = n-i-1;
00314     RealScalar beta;
00315     Scalar h;
00316     matA.col(i).tail(remainingSize).makeHouseholderInPlace(h, beta);
00317     matA.col(i).coeffRef(i+1) = beta;
00318     hCoeffs.coeffRef(i) = h;
00319 
00320     // Apply similarity transformation to remaining columns,
00321     // i.e., compute A = H A H'
00322 
00323     // A = H A
00324     matA.bottomRightCorner(remainingSize, remainingSize)
00325         .applyHouseholderOnTheLeft(matA.col(i).tail(remainingSize-1), h, &temp.coeffRef(0));
00326 
00327     // A = A H'
00328     matA.rightCols(remainingSize)
00329         .applyHouseholderOnTheRight(matA.col(i).tail(remainingSize-1).conjugate(), internal::conj(h), &temp.coeffRef(0));
00330   }
00331 }
00332 
00333 namespace internal {
00334 
00350 template<typename MatrixType> struct HessenbergDecompositionMatrixHReturnType
00351 : public ReturnByValue<HessenbergDecompositionMatrixHReturnType<MatrixType> >
00352 {
00353     typedef typename MatrixType::Index Index;
00354   public:
00359     HessenbergDecompositionMatrixHReturnType(const HessenbergDecomposition<MatrixType>& hess) : m_hess(hess) { }
00360 
00366     template <typename ResultType>
00367     inline void evalTo(ResultType& result) const
00368     {
00369       result = m_hess.packedMatrix();
00370       Index n = result.rows();
00371       if (n>2)
00372         result.bottomLeftCorner(n-2, n-2).template triangularView<Lower>().setZero();
00373     }
00374 
00375     Index rows() const { return m_hess.packedMatrix().rows(); }
00376     Index cols() const { return m_hess.packedMatrix().cols(); }
00377 
00378   protected:
00379     const HessenbergDecomposition<MatrixType>& m_hess;
00380 };
00381 
00382 }
00383 
00384 #endif // EIGEN_HESSENBERGDECOMPOSITION_H


libicr
Author(s): Robert Krug
autogenerated on Mon Jan 6 2014 11:32:46