HouseholderQR.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-2010 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 // Copyright (C) 2010 Vincent Lejeune
00007 //
00008 // This Source Code Form is subject to the terms of the Mozilla
00009 // Public License v. 2.0. If a copy of the MPL was not distributed
00010 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00011 
00012 #ifndef EIGEN_QR_H
00013 #define EIGEN_QR_H
00014 
00015 namespace Eigen { 
00016 
00042 template<typename _MatrixType> class HouseholderQR
00043 {
00044   public:
00045 
00046     typedef _MatrixType MatrixType;
00047     enum {
00048       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00049       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00050       Options = MatrixType::Options,
00051       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00052       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00053     };
00054     typedef typename MatrixType::Scalar Scalar;
00055     typedef typename MatrixType::RealScalar RealScalar;
00056     typedef typename MatrixType::Index Index;
00057     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00058     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00059     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00060     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
00061 
00068     HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
00069 
00076     HouseholderQR(Index rows, Index cols)
00077       : m_qr(rows, cols),
00078         m_hCoeffs((std::min)(rows,cols)),
00079         m_temp(cols),
00080         m_isInitialized(false) {}
00081 
00094     HouseholderQR(const MatrixType& matrix)
00095       : m_qr(matrix.rows(), matrix.cols()),
00096         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00097         m_temp(matrix.cols()),
00098         m_isInitialized(false)
00099     {
00100       compute(matrix);
00101     }
00102 
00120     template<typename Rhs>
00121     inline const internal::solve_retval<HouseholderQR, Rhs>
00122     solve(const MatrixBase<Rhs>& b) const
00123     {
00124       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00125       return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
00126     }
00127 
00136     HouseholderSequenceType householderQ() const
00137     {
00138       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00139       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00140     }
00141 
00145     const MatrixType& matrixQR() const
00146     {
00147         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00148         return m_qr;
00149     }
00150 
00151     HouseholderQR& compute(const MatrixType& matrix);
00152 
00166     typename MatrixType::RealScalar absDeterminant() const;
00167 
00180     typename MatrixType::RealScalar logAbsDeterminant() const;
00181 
00182     inline Index rows() const { return m_qr.rows(); }
00183     inline Index cols() const { return m_qr.cols(); }
00184     
00189     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00190 
00191   protected:
00192     MatrixType m_qr;
00193     HCoeffsType m_hCoeffs;
00194     RowVectorType m_temp;
00195     bool m_isInitialized;
00196 };
00197 
00198 template<typename MatrixType>
00199 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00200 {
00201   using std::abs;
00202   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00203   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00204   return abs(m_qr.diagonal().prod());
00205 }
00206 
00207 template<typename MatrixType>
00208 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00209 {
00210   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00211   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00212   return m_qr.diagonal().cwiseAbs().array().log().sum();
00213 }
00214 
00215 namespace internal {
00216 
00218 template<typename MatrixQR, typename HCoeffs>
00219 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00220 {
00221   typedef typename MatrixQR::Index Index;
00222   typedef typename MatrixQR::Scalar Scalar;
00223   typedef typename MatrixQR::RealScalar RealScalar;
00224   Index rows = mat.rows();
00225   Index cols = mat.cols();
00226   Index size = (std::min)(rows,cols);
00227 
00228   eigen_assert(hCoeffs.size() == size);
00229 
00230   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00231   TempType tempVector;
00232   if(tempData==0)
00233   {
00234     tempVector.resize(cols);
00235     tempData = tempVector.data();
00236   }
00237 
00238   for(Index k = 0; k < size; ++k)
00239   {
00240     Index remainingRows = rows - k;
00241     Index remainingCols = cols - k - 1;
00242 
00243     RealScalar beta;
00244     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00245     mat.coeffRef(k,k) = beta;
00246 
00247     // apply H to remaining part of m_qr from the left
00248     mat.bottomRightCorner(remainingRows, remainingCols)
00249         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00250   }
00251 }
00252 
00254 template<typename MatrixQR, typename HCoeffs>
00255 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
00256                                        typename MatrixQR::Index maxBlockSize=32,
00257                                        typename MatrixQR::Scalar* tempData = 0)
00258 {
00259   typedef typename MatrixQR::Index Index;
00260   typedef typename MatrixQR::Scalar Scalar;
00261   typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00262 
00263   Index rows = mat.rows();
00264   Index cols = mat.cols();
00265   Index size = (std::min)(rows, cols);
00266 
00267   typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00268   TempType tempVector;
00269   if(tempData==0)
00270   {
00271     tempVector.resize(cols);
00272     tempData = tempVector.data();
00273   }
00274 
00275   Index blockSize = (std::min)(maxBlockSize,size);
00276 
00277   Index k = 0;
00278   for (k = 0; k < size; k += blockSize)
00279   {
00280     Index bs = (std::min)(size-k,blockSize);  // actual size of the block
00281     Index tcols = cols - k - bs;            // trailing columns
00282     Index brows = rows-k;                   // rows of the block
00283 
00284     // partition the matrix:
00285     //        A00 | A01 | A02
00286     // mat  = A10 | A11 | A12
00287     //        A20 | A21 | A22
00288     // and performs the qr dec of [A11^T A12^T]^T
00289     // and update [A21^T A22^T]^T using level 3 operations.
00290     // Finally, the algorithm continue on A22
00291 
00292     BlockType A11_21 = mat.block(k,k,brows,bs);
00293     Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00294 
00295     householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00296 
00297     if(tcols)
00298     {
00299       BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00300       apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
00301     }
00302   }
00303 }
00304 
00305 template<typename _MatrixType, typename Rhs>
00306 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
00307   : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
00308 {
00309   EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
00310 
00311   template<typename Dest> void evalTo(Dest& dst) const
00312   {
00313     const Index rows = dec().rows(), cols = dec().cols();
00314     const Index rank = (std::min)(rows, cols);
00315     eigen_assert(rhs().rows() == rows);
00316 
00317     typename Rhs::PlainObject c(rhs());
00318 
00319     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00320     c.applyOnTheLeft(householderSequence(
00321       dec().matrixQR().leftCols(rank),
00322       dec().hCoeffs().head(rank)).transpose()
00323     );
00324 
00325     dec().matrixQR()
00326        .topLeftCorner(rank, rank)
00327        .template triangularView<Upper>()
00328        .solveInPlace(c.topRows(rank));
00329 
00330     dst.topRows(rank) = c.topRows(rank);
00331     dst.bottomRows(cols-rank).setZero();
00332   }
00333 };
00334 
00335 } // end namespace internal
00336 
00343 template<typename MatrixType>
00344 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00345 {
00346   Index rows = matrix.rows();
00347   Index cols = matrix.cols();
00348   Index size = (std::min)(rows,cols);
00349 
00350   m_qr = matrix;
00351   m_hCoeffs.resize(size);
00352 
00353   m_temp.resize(cols);
00354 
00355   internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
00356 
00357   m_isInitialized = true;
00358   return *this;
00359 }
00360 
00365 template<typename Derived>
00366 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00367 MatrixBase<Derived>::householderQr() const
00368 {
00369   return HouseholderQR<PlainObject>(eval());
00370 }
00371 
00372 } // end namespace Eigen
00373 
00374 #endif // EIGEN_QR_H


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Sat Jun 8 2019 19:37:24