Eigen  3.3.3
HouseholderQR.h
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 
00044 template<typename _MatrixType> class HouseholderQR
00045 {
00046   public:
00047 
00048     typedef _MatrixType MatrixType;
00049     enum {
00050       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00051       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00052       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00053       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00054     };
00055     typedef typename MatrixType::Scalar Scalar;
00056     typedef typename MatrixType::RealScalar RealScalar;
00057     // FIXME should be int
00058     typedef typename MatrixType::StorageIndex StorageIndex;
00059     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, (MatrixType::Flags&RowMajorBit) ? RowMajor : ColMajor, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00060     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00061     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00062     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
00063 
00070     HouseholderQR() : m_qr(), m_hCoeffs(), m_temp(), m_isInitialized(false) {}
00071 
00078     HouseholderQR(Index rows, Index cols)
00079       : m_qr(rows, cols),
00080         m_hCoeffs((std::min)(rows,cols)),
00081         m_temp(cols),
00082         m_isInitialized(false) {}
00083 
00096     template<typename InputType>
00097     explicit HouseholderQR(const EigenBase<InputType>& matrix)
00098       : m_qr(matrix.rows(), matrix.cols()),
00099         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00100         m_temp(matrix.cols()),
00101         m_isInitialized(false)
00102     {
00103       compute(matrix.derived());
00104     }
00105 
00106 
00114     template<typename InputType>
00115     explicit HouseholderQR(EigenBase<InputType>& matrix)
00116       : m_qr(matrix.derived()),
00117         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00118         m_temp(matrix.cols()),
00119         m_isInitialized(false)
00120     {
00121       computeInPlace();
00122     }
00123 
00138     template<typename Rhs>
00139     inline const Solve<HouseholderQR, Rhs>
00140     solve(const MatrixBase<Rhs>& b) const
00141     {
00142       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00143       return Solve<HouseholderQR, Rhs>(*this, b.derived());
00144     }
00145 
00154     HouseholderSequenceType householderQ() const
00155     {
00156       eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00157       return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
00158     }
00159 
00163     const MatrixType& matrixQR() const
00164     {
00165         eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00166         return m_qr;
00167     }
00168 
00169     template<typename InputType>
00170     HouseholderQR& compute(const EigenBase<InputType>& matrix) {
00171       m_qr = matrix.derived();
00172       computeInPlace();
00173       return *this;
00174     }
00175 
00189     typename MatrixType::RealScalar absDeterminant() const;
00190 
00203     typename MatrixType::RealScalar logAbsDeterminant() const;
00204 
00205     inline Index rows() const { return m_qr.rows(); }
00206     inline Index cols() const { return m_qr.cols(); }
00207     
00212     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
00213     
00214     #ifndef EIGEN_PARSED_BY_DOXYGEN
00215     template<typename RhsType, typename DstType>
00216     EIGEN_DEVICE_FUNC
00217     void _solve_impl(const RhsType &rhs, DstType &dst) const;
00218     #endif
00219 
00220   protected:
00221     
00222     static void check_template_parameters()
00223     {
00224       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00225     }
00226 
00227     void computeInPlace();
00228     
00229     MatrixType m_qr;
00230     HCoeffsType m_hCoeffs;
00231     RowVectorType m_temp;
00232     bool m_isInitialized;
00233 };
00234 
00235 template<typename MatrixType>
00236 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const
00237 {
00238   using std::abs;
00239   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00240   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00241   return abs(m_qr.diagonal().prod());
00242 }
00243 
00244 template<typename MatrixType>
00245 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const
00246 {
00247   eigen_assert(m_isInitialized && "HouseholderQR is not initialized.");
00248   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00249   return m_qr.diagonal().cwiseAbs().array().log().sum();
00250 }
00251 
00252 namespace internal {
00253 
00255 template<typename MatrixQR, typename HCoeffs>
00256 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0)
00257 {
00258   typedef typename MatrixQR::Scalar Scalar;
00259   typedef typename MatrixQR::RealScalar RealScalar;
00260   Index rows = mat.rows();
00261   Index cols = mat.cols();
00262   Index size = (std::min)(rows,cols);
00263 
00264   eigen_assert(hCoeffs.size() == size);
00265 
00266   typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType;
00267   TempType tempVector;
00268   if(tempData==0)
00269   {
00270     tempVector.resize(cols);
00271     tempData = tempVector.data();
00272   }
00273 
00274   for(Index k = 0; k < size; ++k)
00275   {
00276     Index remainingRows = rows - k;
00277     Index remainingCols = cols - k - 1;
00278 
00279     RealScalar beta;
00280     mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta);
00281     mat.coeffRef(k,k) = beta;
00282 
00283     // apply H to remaining part of m_qr from the left
00284     mat.bottomRightCorner(remainingRows, remainingCols)
00285         .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1);
00286   }
00287 }
00288 
00290 template<typename MatrixQR, typename HCoeffs,
00291   typename MatrixQRScalar = typename MatrixQR::Scalar,
00292   bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)>
00293 struct householder_qr_inplace_blocked
00294 {
00295   // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h
00296   static void run(MatrixQR& mat, HCoeffs& hCoeffs, Index maxBlockSize=32,
00297       typename MatrixQR::Scalar* tempData = 0)
00298   {
00299     typedef typename MatrixQR::Scalar Scalar;
00300     typedef Block<MatrixQR,Dynamic,Dynamic> BlockType;
00301 
00302     Index rows = mat.rows();
00303     Index cols = mat.cols();
00304     Index size = (std::min)(rows, cols);
00305 
00306     typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType;
00307     TempType tempVector;
00308     if(tempData==0)
00309     {
00310       tempVector.resize(cols);
00311       tempData = tempVector.data();
00312     }
00313 
00314     Index blockSize = (std::min)(maxBlockSize,size);
00315 
00316     Index k = 0;
00317     for (k = 0; k < size; k += blockSize)
00318     {
00319       Index bs = (std::min)(size-k,blockSize);  // actual size of the block
00320       Index tcols = cols - k - bs;              // trailing columns
00321       Index brows = rows-k;                     // rows of the block
00322 
00323       // partition the matrix:
00324       //        A00 | A01 | A02
00325       // mat  = A10 | A11 | A12
00326       //        A20 | A21 | A22
00327       // and performs the qr dec of [A11^T A12^T]^T
00328       // and update [A21^T A22^T]^T using level 3 operations.
00329       // Finally, the algorithm continue on A22
00330 
00331       BlockType A11_21 = mat.block(k,k,brows,bs);
00332       Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs);
00333 
00334       householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData);
00335 
00336       if(tcols)
00337       {
00338         BlockType A21_22 = mat.block(k,k+bs,brows,tcols);
00339         apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment, false); // false == backward
00340       }
00341     }
00342   }
00343 };
00344 
00345 } // end namespace internal
00346 
00347 #ifndef EIGEN_PARSED_BY_DOXYGEN
00348 template<typename _MatrixType>
00349 template<typename RhsType, typename DstType>
00350 void HouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
00351 {
00352   const Index rank = (std::min)(rows(), cols());
00353   eigen_assert(rhs.rows() == rows());
00354 
00355   typename RhsType::PlainObject c(rhs);
00356 
00357   // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00358   c.applyOnTheLeft(householderSequence(
00359     m_qr.leftCols(rank),
00360     m_hCoeffs.head(rank)).transpose()
00361   );
00362 
00363   m_qr.topLeftCorner(rank, rank)
00364       .template triangularView<Upper>()
00365       .solveInPlace(c.topRows(rank));
00366 
00367   dst.topRows(rank) = c.topRows(rank);
00368   dst.bottomRows(cols()-rank).setZero();
00369 }
00370 #endif
00371 
00378 template<typename MatrixType>
00379 void HouseholderQR<MatrixType>::computeInPlace()
00380 {
00381   check_template_parameters();
00382   
00383   Index rows = m_qr.rows();
00384   Index cols = m_qr.cols();
00385   Index size = (std::min)(rows,cols);
00386 
00387   m_hCoeffs.resize(size);
00388 
00389   m_temp.resize(cols);
00390 
00391   internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data());
00392 
00393   m_isInitialized = true;
00394 }
00395 
00400 template<typename Derived>
00401 const HouseholderQR<typename MatrixBase<Derived>::PlainObject>
00402 MatrixBase<Derived>::householderQr() const
00403 {
00404   return HouseholderQR<PlainObject>(eval());
00405 }
00406 
00407 } // end namespace Eigen
00408 
00409 #endif // EIGEN_QR_H
 All Classes Functions Variables Typedefs Enumerations Enumerator Friends