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 typename HouseholderSequence<MatrixType,HCoeffsType>::ConjugateReturnType 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 
00082     HouseholderQR(const MatrixType& matrix)
00083       : m_qr(matrix.rows(), matrix.cols()),
00084         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00085         m_temp(matrix.cols()),
00086         m_isInitialized(false)
00087     {
00088       compute(matrix);
00089     }
00090 
00108     template<typename Rhs>
00109     inline const internal::solve_retval<HouseholderQR, Rhs>
00110     solve(const MatrixBase<Rhs>& b) const
00111     {
00112       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00113       return internal::solve_retval<HouseholderQR, Rhs>(*this, b.derived());
00114     }
00115 
00116     HouseholderSequenceType householderQ() const
00117     {
00118       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00119       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00120     }
00121 
00125     const MatrixType& matrixQR() const
00126     {
00127         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00128         return m_qr;
00129     }
00130 
00131     HouseholderQR& compute(const MatrixType& matrix);
00132 
00146     typename MatrixType::RealScalar absDeterminant() const;
00147 
00160     typename MatrixType::RealScalar logAbsDeterminant() const;
00161 
00162     inline Index rows() const { return m_qr.rows(); }
00163     inline Index cols() const { return m_qr.cols(); }
00164     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00165 
00166   protected:
00167     MatrixType m_qr;
00168     HCoeffsType m_hCoeffs;
00169     RowVectorType m_temp;
00170     bool m_isInitialized;
00171 };
00172 
00173 template<typename MatrixType>
00174 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00175 {
00176   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00177   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00178   return internal::abs(m_qr.diagonal().prod());
00179 }
00180 
00181 template<typename MatrixType>
00182 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00183 {
00184   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00185   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00186   return m_qr.diagonal().cwiseAbs().array().log().sum();
00187 }
00188 
00189 namespace internal {
00190 
00192 template<typename MatrixQR, typename HCoeffs>
00193 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00194 {
00195   typedef typename MatrixQR::Index Index;
00196   typedef typename MatrixQR::Scalar Scalar;
00197   typedef typename MatrixQR::RealScalar RealScalar;
00198   Index rows = mat.rows();
00199   Index cols = mat.cols();
00200   Index size = (std::min)(rows,cols);
00201 
00202   eigen_assert(hCoeffs.size() == size);
00203 
00204   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00205   TempType tempVector;
00206   if(tempData==0)
00207   {
00208     tempVector.resize(cols);
00209     tempData = tempVector.data();
00210   }
00211 
00212   for(Index k = 0; k < size; ++k)
00213   {
00214     Index remainingRows = rows - k;
00215     Index remainingCols = cols - k - 1;
00216 
00217     RealScalar beta;
00218     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00219     mat.coeffRef(k,k) = beta;
00220 
00221     // apply H to remaining part of m_qr from the left
00222     mat.bottomRightCorner(remainingRows, remainingCols)
00223         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00224   }
00225 }
00226 
00228 template<typename MatrixQR, typename HCoeffs>
00229 void householder_qr_inplace_blocked(MatrixQR& mat, HCoeffs& hCoeffs,
00230                                        typename MatrixQR::Index maxBlockSize=32,
00231                                        typename MatrixQR::Scalar* tempData = 0)
00232 {
00233   typedef typename MatrixQR::Index Index;
00234   typedef typename MatrixQR::Scalar Scalar;
00235   typedef typename MatrixQR::RealScalar RealScalar;
00236   typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00237 
00238   Index rows = mat.rows();
00239   Index cols = mat.cols();
00240   Index size = (std::min)(rows, cols);
00241 
00242   typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00243   TempType tempVector;
00244   if(tempData==0)
00245   {
00246     tempVector.resize(cols);
00247     tempData = tempVector.data();
00248   }
00249 
00250   Index blockSize = (std::min)(maxBlockSize,size);
00251 
00252   Index k = 0;
00253   for (k = 0; k < size; k += blockSize)
00254   {
00255     Index bs = (std::min)(size-k,blockSize);  // actual size of the block
00256     Index tcols = cols - k - bs;            // trailing columns
00257     Index brows = rows-k;                   // rows of the block
00258 
00259     // partition the matrix:
00260     //        A00 | A01 | A02
00261     // mat  = A10 | A11 | A12
00262     //        A20 | A21 | A22
00263     // and performs the qr dec of [A11^T A12^T]^T
00264     // and update [A21^T A22^T]^T using level 3 operations.
00265     // Finally, the algorithm continue on A22
00266 
00267     BlockType A11_21 = mat.block(k,k,brows,bs);
00268     Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00269 
00270     householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00271 
00272     if(tcols)
00273     {
00274       BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00275       apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint());
00276     }
00277   }
00278 }
00279 
00280 template<typename _MatrixType, typename Rhs>
00281 struct solve_retval<HouseholderQR<_MatrixType>, Rhs>
00282   : solve_retval_base<HouseholderQR<_MatrixType>, Rhs>
00283 {
00284   EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs)
00285 
00286   template<typename Dest> void evalTo(Dest& dst) const
00287   {
00288     const Index rows = dec().rows(), cols = dec().cols();
00289     const Index rank = (std::min)(rows, cols);
00290     eigen_assert(rhs().rows() == rows);
00291 
00292     typename Rhs::PlainObject c(rhs());
00293 
00294     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00295     c.applyOnTheLeft(householderSequence(
00296       dec().matrixQR().leftCols(rank),
00297       dec().hCoeffs().head(rank)).transpose()
00298     );
00299 
00300     dec().matrixQR()
00301        .topLeftCorner(rank, rank)
00302        .template triangularView<Upper>()
00303        .solveInPlace(c.topRows(rank));
00304 
00305     dst.topRows(rank) = c.topRows(rank);
00306     dst.bottomRows(cols-rank).setZero();
00307   }
00308 };
00309 
00310 } // end namespace internal
00311 
00312 template<typename MatrixType>
00313 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00314 {
00315   Index rows = matrix.rows();
00316   Index cols = matrix.cols();
00317   Index size = (std::min)(rows,cols);
00318 
00319   m_qr = matrix;
00320   m_hCoeffs.resize(size);
00321 
00322   m_temp.resize(cols);
00323 
00324   internal::householder_qr_inplace_blocked(m_qr, m_hCoeffs, 48, m_temp.data());
00325 
00326   m_isInitialized = true;
00327   return *this;
00328 }
00329 
00334 template<typename Derived>
00335 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00336 MatrixBase<Derived>::householderQr() const
00337 {
00338   return HouseholderQR<PlainObject>(eval());
00339 }
00340 
00341 } // end namespace Eigen
00342 
00343 #endif // EIGEN_QR_H


win_eigen
Author(s): Daniel Stonier
autogenerated on Mon Oct 6 2014 12:24:42