$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
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 00193 static void check_template_parameters() 00194 { 00195 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 00196 } 00197 00198 MatrixType m_qr; 00199 HCoeffsType m_hCoeffs; 00200 RowVectorType m_temp; 00201 bool m_isInitialized; 00202 }; 00203 00204 template<typename MatrixType> 00205 typename MatrixType::RealScalar HouseholderQR<MatrixType>::absDeterminant() const 00206 { 00207 using std::abs; 00208 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 00209 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00210 return abs(m_qr.diagonal().prod()); 00211 } 00212 00213 template<typename MatrixType> 00214 typename MatrixType::RealScalar HouseholderQR<MatrixType>::logAbsDeterminant() const 00215 { 00216 eigen_assert(m_isInitialized && "HouseholderQR is not initialized."); 00217 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 00218 return m_qr.diagonal().cwiseAbs().array().log().sum(); 00219 } 00220 00221 namespace internal { 00222 00224 template<typename MatrixQR, typename HCoeffs> 00225 void householder_qr_inplace_unblocked(MatrixQR& mat, HCoeffs& hCoeffs, typename MatrixQR::Scalar* tempData = 0) 00226 { 00227 typedef typename MatrixQR::Index Index; 00228 typedef typename MatrixQR::Scalar Scalar; 00229 typedef typename MatrixQR::RealScalar RealScalar; 00230 Index rows = mat.rows(); 00231 Index cols = mat.cols(); 00232 Index size = (std::min)(rows,cols); 00233 00234 eigen_assert(hCoeffs.size() == size); 00235 00236 typedef Matrix<Scalar,MatrixQR::ColsAtCompileTime,1> TempType; 00237 TempType tempVector; 00238 if(tempData==0) 00239 { 00240 tempVector.resize(cols); 00241 tempData = tempVector.data(); 00242 } 00243 00244 for(Index k = 0; k < size; ++k) 00245 { 00246 Index remainingRows = rows - k; 00247 Index remainingCols = cols - k - 1; 00248 00249 RealScalar beta; 00250 mat.col(k).tail(remainingRows).makeHouseholderInPlace(hCoeffs.coeffRef(k), beta); 00251 mat.coeffRef(k,k) = beta; 00252 00253 // apply H to remaining part of m_qr from the left 00254 mat.bottomRightCorner(remainingRows, remainingCols) 00255 .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), hCoeffs.coeffRef(k), tempData+k+1); 00256 } 00257 } 00258 00260 template<typename MatrixQR, typename HCoeffs, 00261 typename MatrixQRScalar = typename MatrixQR::Scalar, 00262 bool InnerStrideIsOne = (MatrixQR::InnerStrideAtCompileTime == 1 && HCoeffs::InnerStrideAtCompileTime == 1)> 00263 struct householder_qr_inplace_blocked 00264 { 00265 // This is specialized for MKL-supported Scalar types in HouseholderQR_MKL.h 00266 static void run(MatrixQR& mat, HCoeffs& hCoeffs, 00267 typename MatrixQR::Index maxBlockSize=32, 00268 typename MatrixQR::Scalar* tempData = 0) 00269 { 00270 typedef typename MatrixQR::Index Index; 00271 typedef typename MatrixQR::Scalar Scalar; 00272 typedef Block<MatrixQR,Dynamic,Dynamic> BlockType; 00273 00274 Index rows = mat.rows(); 00275 Index cols = mat.cols(); 00276 Index size = (std::min)(rows, cols); 00277 00278 typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixQR::MaxColsAtCompileTime,1> TempType; 00279 TempType tempVector; 00280 if(tempData==0) 00281 { 00282 tempVector.resize(cols); 00283 tempData = tempVector.data(); 00284 } 00285 00286 Index blockSize = (std::min)(maxBlockSize,size); 00287 00288 Index k = 0; 00289 for (k = 0; k < size; k += blockSize) 00290 { 00291 Index bs = (std::min)(size-k,blockSize); // actual size of the block 00292 Index tcols = cols - k - bs; // trailing columns 00293 Index brows = rows-k; // rows of the block 00294 00295 // partition the matrix: 00296 // A00 | A01 | A02 00297 // mat = A10 | A11 | A12 00298 // A20 | A21 | A22 00299 // and performs the qr dec of [A11^T A12^T]^T 00300 // and update [A21^T A22^T]^T using level 3 operations. 00301 // Finally, the algorithm continue on A22 00302 00303 BlockType A11_21 = mat.block(k,k,brows,bs); 00304 Block<HCoeffs,Dynamic,1> hCoeffsSegment = hCoeffs.segment(k,bs); 00305 00306 householder_qr_inplace_unblocked(A11_21, hCoeffsSegment, tempData); 00307 00308 if(tcols) 00309 { 00310 BlockType A21_22 = mat.block(k,k+bs,brows,tcols); 00311 apply_block_householder_on_the_left(A21_22,A11_21,hCoeffsSegment.adjoint()); 00312 } 00313 } 00314 } 00315 }; 00316 00317 template<typename _MatrixType, typename Rhs> 00318 struct solve_retval<HouseholderQR<_MatrixType>, Rhs> 00319 : solve_retval_base<HouseholderQR<_MatrixType>, Rhs> 00320 { 00321 EIGEN_MAKE_SOLVE_HELPERS(HouseholderQR<_MatrixType>,Rhs) 00322 00323 template<typename Dest> void evalTo(Dest& dst) const 00324 { 00325 const Index rows = dec().rows(), cols = dec().cols(); 00326 const Index rank = (std::min)(rows, cols); 00327 eigen_assert(rhs().rows() == rows); 00328 00329 typename Rhs::PlainObject c(rhs()); 00330 00331 // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T 00332 c.applyOnTheLeft(householderSequence( 00333 dec().matrixQR().leftCols(rank), 00334 dec().hCoeffs().head(rank)).transpose() 00335 ); 00336 00337 dec().matrixQR() 00338 .topLeftCorner(rank, rank) 00339 .template triangularView<Upper>() 00340 .solveInPlace(c.topRows(rank)); 00341 00342 dst.topRows(rank) = c.topRows(rank); 00343 dst.bottomRows(cols-rank).setZero(); 00344 } 00345 }; 00346 00347 } // end namespace internal 00348 00355 template<typename MatrixType> 00356 HouseholderQR<MatrixType>& HouseholderQR<MatrixType>::compute(const MatrixType& matrix) 00357 { 00358 check_template_parameters(); 00359 00360 Index rows = matrix.rows(); 00361 Index cols = matrix.cols(); 00362 Index size = (std::min)(rows,cols); 00363 00364 m_qr = matrix; 00365 m_hCoeffs.resize(size); 00366 00367 m_temp.resize(cols); 00368 00369 internal::householder_qr_inplace_blocked<MatrixType, HCoeffsType>::run(m_qr, m_hCoeffs, 48, m_temp.data()); 00370 00371 m_isInitialized = true; 00372 return *this; 00373 } 00374 00379 template<typename Derived> 00380 const HouseholderQR<typename MatrixBase<Derived>::PlainObject> 00381 MatrixBase<Derived>::householderQr() const 00382 { 00383 return HouseholderQR<PlainObject>(eval()); 00384 } 00385 00386 } // end namespace Eigen 00387 00388 #endif // EIGEN_QR_H