![]() |
Eigen
3.3.3
|
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