$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) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam@inria.fr> 00005 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 00010 00011 #ifndef EIGEN_SPARSE_QR_H 00012 #define EIGEN_SPARSE_QR_H 00013 00014 namespace Eigen { 00015 00016 template<typename MatrixType, typename OrderingType> class SparseQR; 00017 template<typename SparseQRType> struct SparseQRMatrixQReturnType; 00018 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType; 00019 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct; 00020 namespace internal { 00021 template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> > 00022 { 00023 typedef typename SparseQRType::MatrixType ReturnType; 00024 typedef typename ReturnType::Index Index; 00025 typedef typename ReturnType::StorageKind StorageKind; 00026 }; 00027 template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> > 00028 { 00029 typedef typename SparseQRType::MatrixType ReturnType; 00030 }; 00031 template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> > 00032 { 00033 typedef typename Derived::PlainObject ReturnType; 00034 }; 00035 } // End namespace internal 00036 00064 template<typename _MatrixType, typename _OrderingType> 00065 class SparseQR 00066 { 00067 public: 00068 typedef _MatrixType MatrixType; 00069 typedef _OrderingType OrderingType; 00070 typedef typename MatrixType::Scalar Scalar; 00071 typedef typename MatrixType::RealScalar RealScalar; 00072 typedef typename MatrixType::Index Index; 00073 typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType; 00074 typedef Matrix<Index, Dynamic, 1> IndexVector; 00075 typedef Matrix<Scalar, Dynamic, 1> ScalarVector; 00076 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 00077 public: 00078 SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) 00079 { } 00080 00087 SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) 00088 { 00089 compute(mat); 00090 } 00091 00098 void compute(const MatrixType& mat) 00099 { 00100 analyzePattern(mat); 00101 factorize(mat); 00102 } 00103 void analyzePattern(const MatrixType& mat); 00104 void factorize(const MatrixType& mat); 00105 00108 inline Index rows() const { return m_pmat.rows(); } 00109 00112 inline Index cols() const { return m_pmat.cols();} 00113 00116 const QRMatrixType& matrixR() const { return m_R; } 00117 00122 Index rank() const 00123 { 00124 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00125 return m_nonzeropivots; 00126 } 00127 00146 SparseQRMatrixQReturnType<SparseQR> matrixQ() const 00147 { return SparseQRMatrixQReturnType<SparseQR>(*this); } 00148 00152 const PermutationType& colsPermutation() const 00153 { 00154 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00155 return m_outputPerm_c; 00156 } 00157 00161 std::string lastErrorMessage() const { return m_lastError; } 00162 00164 template<typename Rhs, typename Dest> 00165 bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const 00166 { 00167 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00168 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00169 00170 Index rank = this->rank(); 00171 00172 // Compute Q^T * b; 00173 typename Dest::PlainObject y, b; 00174 y = this->matrixQ().transpose() * B; 00175 b = y; 00176 00177 // Solve with the triangular matrix R 00178 y.resize((std::max)(cols(),Index(y.rows())),y.cols()); 00179 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank)); 00180 y.bottomRows(y.rows()-rank).setZero(); 00181 00182 // Apply the column permutation 00183 if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); 00184 else dest = y.topRows(cols()); 00185 00186 m_info = Success; 00187 return true; 00188 } 00189 00190 00196 void setPivotThreshold(const RealScalar& threshold) 00197 { 00198 m_useDefaultThreshold = false; 00199 m_threshold = threshold; 00200 } 00201 00206 template<typename Rhs> 00207 inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 00208 { 00209 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00210 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00211 return internal::solve_retval<SparseQR, Rhs>(*this, B.derived()); 00212 } 00213 template<typename Rhs> 00214 inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 00215 { 00216 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00217 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00218 return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived()); 00219 } 00220 00229 ComputationInfo info() const 00230 { 00231 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00232 return m_info; 00233 } 00234 00235 protected: 00236 inline void sort_matrix_Q() 00237 { 00238 if(this->m_isQSorted) return; 00239 // The matrix Q is sorted during the transposition 00240 SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q); 00241 this->m_Q = mQrm; 00242 this->m_isQSorted = true; 00243 } 00244 00245 00246 protected: 00247 bool m_isInitialized; 00248 bool m_analysisIsok; 00249 bool m_factorizationIsok; 00250 mutable ComputationInfo m_info; 00251 std::string m_lastError; 00252 QRMatrixType m_pmat; // Temporary matrix 00253 QRMatrixType m_R; // The triangular factor matrix 00254 QRMatrixType m_Q; // The orthogonal reflectors 00255 ScalarVector m_hcoeffs; // The Householder coefficients 00256 PermutationType m_perm_c; // Fill-reducing Column permutation 00257 PermutationType m_pivotperm; // The permutation for rank revealing 00258 PermutationType m_outputPerm_c; // The final column permutation 00259 RealScalar m_threshold; // Threshold to determine null Householder reflections 00260 bool m_useDefaultThreshold; // Use default threshold 00261 Index m_nonzeropivots; // Number of non zero pivots found 00262 IndexVector m_etree; // Column elimination tree 00263 IndexVector m_firstRowElt; // First element in each row 00264 bool m_isQSorted; // whether Q is sorted or not 00265 bool m_isEtreeOk; // whether the elimination tree match the initial input matrix 00266 00267 template <typename, typename > friend struct SparseQR_QProduct; 00268 template <typename > friend struct SparseQRMatrixQReturnType; 00269 00270 }; 00271 00281 template <typename MatrixType, typename OrderingType> 00282 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) 00283 { 00284 eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); 00285 // Copy to a column major matrix if the input is rowmajor 00286 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat); 00287 // Compute the column fill reducing ordering 00288 OrderingType ord; 00289 ord(matCpy, m_perm_c); 00290 Index n = mat.cols(); 00291 Index m = mat.rows(); 00292 Index diagSize = (std::min)(m,n); 00293 00294 if (!m_perm_c.size()) 00295 { 00296 m_perm_c.resize(n); 00297 m_perm_c.indices().setLinSpaced(n, 0,n-1); 00298 } 00299 00300 // Compute the column elimination tree of the permuted matrix 00301 m_outputPerm_c = m_perm_c.inverse(); 00302 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 00303 m_isEtreeOk = true; 00304 00305 m_R.resize(m, n); 00306 m_Q.resize(m, diagSize); 00307 00308 // Allocate space for nonzero elements : rough estimation 00309 m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree 00310 m_Q.reserve(2*mat.nonZeros()); 00311 m_hcoeffs.resize(diagSize); 00312 m_analysisIsok = true; 00313 } 00314 00322 template <typename MatrixType, typename OrderingType> 00323 void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) 00324 { 00325 using std::abs; 00326 using std::max; 00327 00328 eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); 00329 Index m = mat.rows(); 00330 Index n = mat.cols(); 00331 Index diagSize = (std::min)(m,n); 00332 IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes 00333 IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q 00334 Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q 00335 ScalarVector tval(m); // The dense vector used to compute the current column 00336 RealScalar pivotThreshold = m_threshold; 00337 00338 m_R.setZero(); 00339 m_Q.setZero(); 00340 m_pmat = mat; 00341 if(!m_isEtreeOk) 00342 { 00343 m_outputPerm_c = m_perm_c.inverse(); 00344 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 00345 m_isEtreeOk = true; 00346 } 00347 00348 m_pmat.uncompress(); // To have the innerNonZeroPtr allocated 00349 00350 // Apply the fill-in reducing permutation lazily: 00351 { 00352 // If the input is row major, copy the original column indices, 00353 // otherwise directly use the input matrix 00354 // 00355 IndexVector originalOuterIndicesCpy; 00356 const Index *originalOuterIndices = mat.outerIndexPtr(); 00357 if(MatrixType::IsRowMajor) 00358 { 00359 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1); 00360 originalOuterIndices = originalOuterIndicesCpy.data(); 00361 } 00362 00363 for (int i = 0; i < n; i++) 00364 { 00365 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; 00366 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; 00367 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; 00368 } 00369 } 00370 00371 /* Compute the default threshold as in MatLab, see: 00372 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing 00373 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 00374 */ 00375 if(m_useDefaultThreshold) 00376 { 00377 RealScalar max2Norm = 0.0; 00378 for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); 00379 if(max2Norm==RealScalar(0)) 00380 max2Norm = RealScalar(1); 00381 pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); 00382 } 00383 00384 // Initialize the numerical permutation 00385 m_pivotperm.setIdentity(n); 00386 00387 Index nonzeroCol = 0; // Record the number of valid pivots 00388 m_Q.startVec(0); 00389 00390 // Left looking rank-revealing QR factorization: compute a column of R and Q at a time 00391 for (Index col = 0; col < n; ++col) 00392 { 00393 mark.setConstant(-1); 00394 m_R.startVec(col); 00395 mark(nonzeroCol) = col; 00396 Qidx(0) = nonzeroCol; 00397 nzcolR = 0; nzcolQ = 1; 00398 bool found_diag = nonzeroCol>=m; 00399 tval.setZero(); 00400 00401 // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., 00402 // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. 00403 // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, 00404 // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. 00405 for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) 00406 { 00407 Index curIdx = nonzeroCol; 00408 if(itp) curIdx = itp.row(); 00409 if(curIdx == nonzeroCol) found_diag = true; 00410 00411 // Get the nonzeros indexes of the current column of R 00412 Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here 00413 if (st < 0 ) 00414 { 00415 m_lastError = "Empty row found during numerical factorization"; 00416 m_info = InvalidInput; 00417 return; 00418 } 00419 00420 // Traverse the etree 00421 Index bi = nzcolR; 00422 for (; mark(st) != col; st = m_etree(st)) 00423 { 00424 Ridx(nzcolR) = st; // Add this row to the list, 00425 mark(st) = col; // and mark this row as visited 00426 nzcolR++; 00427 } 00428 00429 // Reverse the list to get the topological ordering 00430 Index nt = nzcolR-bi; 00431 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1)); 00432 00433 // Copy the current (curIdx,pcol) value of the input matrix 00434 if(itp) tval(curIdx) = itp.value(); 00435 else tval(curIdx) = Scalar(0); 00436 00437 // Compute the pattern of Q(:,k) 00438 if(curIdx > nonzeroCol && mark(curIdx) != col ) 00439 { 00440 Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q, 00441 mark(curIdx) = col; // and mark it as visited 00442 nzcolQ++; 00443 } 00444 } 00445 00446 // Browse all the indexes of R(:,col) in reverse order 00447 for (Index i = nzcolR-1; i >= 0; i--) 00448 { 00449 Index curIdx = Ridx(i); 00450 00451 // Apply the curIdx-th householder vector to the current column (temporarily stored into tval) 00452 Scalar tdot(0); 00453 00454 // First compute q' * tval 00455 tdot = m_Q.col(curIdx).dot(tval); 00456 00457 tdot *= m_hcoeffs(curIdx); 00458 00459 // Then update tval = tval - q * tau 00460 // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse") 00461 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 00462 tval(itq.row()) -= itq.value() * tdot; 00463 00464 // Detect fill-in for the current column of Q 00465 if(m_etree(Ridx(i)) == nonzeroCol) 00466 { 00467 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 00468 { 00469 Index iQ = itq.row(); 00470 if (mark(iQ) != col) 00471 { 00472 Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q, 00473 mark(iQ) = col; // and mark it as visited 00474 } 00475 } 00476 } 00477 } // End update current column 00478 00479 Scalar tau = 0; 00480 RealScalar beta = 0; 00481 00482 if(nonzeroCol < diagSize) 00483 { 00484 // Compute the Householder reflection that eliminate the current column 00485 // FIXME this step should call the Householder module. 00486 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); 00487 00488 // First, the squared norm of Q((col+1):m, col) 00489 RealScalar sqrNorm = 0.; 00490 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); 00491 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) 00492 { 00493 beta = numext::real(c0); 00494 tval(Qidx(0)) = 1; 00495 } 00496 else 00497 { 00498 using std::sqrt; 00499 beta = sqrt(numext::abs2(c0) + sqrNorm); 00500 if(numext::real(c0) >= RealScalar(0)) 00501 beta = -beta; 00502 tval(Qidx(0)) = 1; 00503 for (Index itq = 1; itq < nzcolQ; ++itq) 00504 tval(Qidx(itq)) /= (c0 - beta); 00505 tau = numext::conj((beta-c0) / beta); 00506 00507 } 00508 } 00509 00510 // Insert values in R 00511 for (Index i = nzcolR-1; i >= 0; i--) 00512 { 00513 Index curIdx = Ridx(i); 00514 if(curIdx < nonzeroCol) 00515 { 00516 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); 00517 tval(curIdx) = Scalar(0.); 00518 } 00519 } 00520 00521 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold) 00522 { 00523 m_R.insertBackByOuterInner(col, nonzeroCol) = beta; 00524 // The householder coefficient 00525 m_hcoeffs(nonzeroCol) = tau; 00526 // Record the householder reflections 00527 for (Index itq = 0; itq < nzcolQ; ++itq) 00528 { 00529 Index iQ = Qidx(itq); 00530 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ); 00531 tval(iQ) = Scalar(0.); 00532 } 00533 nonzeroCol++; 00534 if(nonzeroCol<diagSize) 00535 m_Q.startVec(nonzeroCol); 00536 } 00537 else 00538 { 00539 // Zero pivot found: move implicitly this column to the end 00540 for (Index j = nonzeroCol; j < n-1; j++) 00541 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]); 00542 00543 // Recompute the column elimination tree 00544 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); 00545 m_isEtreeOk = false; 00546 } 00547 } 00548 00549 m_hcoeffs.tail(diagSize-nonzeroCol).setZero(); 00550 00551 // Finalize the column pointers of the sparse matrices R and Q 00552 m_Q.finalize(); 00553 m_Q.makeCompressed(); 00554 m_R.finalize(); 00555 m_R.makeCompressed(); 00556 m_isQSorted = false; 00557 00558 m_nonzeropivots = nonzeroCol; 00559 00560 if(nonzeroCol<n) 00561 { 00562 // Permute the triangular factor to put the 'dead' columns to the end 00563 QRMatrixType tempR(m_R); 00564 m_R = tempR * m_pivotperm; 00565 00566 // Update the column permutation 00567 m_outputPerm_c = m_outputPerm_c * m_pivotperm; 00568 } 00569 00570 m_isInitialized = true; 00571 m_factorizationIsok = true; 00572 m_info = Success; 00573 } 00574 00575 namespace internal { 00576 00577 template<typename _MatrixType, typename OrderingType, typename Rhs> 00578 struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs> 00579 : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs> 00580 { 00581 typedef SparseQR<_MatrixType,OrderingType> Dec; 00582 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00583 00584 template<typename Dest> void evalTo(Dest& dst) const 00585 { 00586 dec()._solve(rhs(),dst); 00587 } 00588 }; 00589 template<typename _MatrixType, typename OrderingType, typename Rhs> 00590 struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs> 00591 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs> 00592 { 00593 typedef SparseQR<_MatrixType, OrderingType> Dec; 00594 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs) 00595 00596 template<typename Dest> void evalTo(Dest& dst) const 00597 { 00598 this->defaultEvalTo(dst); 00599 } 00600 }; 00601 } // end namespace internal 00602 00603 template <typename SparseQRType, typename Derived> 00604 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> > 00605 { 00606 typedef typename SparseQRType::QRMatrixType MatrixType; 00607 typedef typename SparseQRType::Scalar Scalar; 00608 typedef typename SparseQRType::Index Index; 00609 // Get the references 00610 SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 00611 m_qr(qr),m_other(other),m_transpose(transpose) {} 00612 inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); } 00613 inline Index cols() const { return m_other.cols(); } 00614 00615 // Assign to a vector 00616 template<typename DesType> 00617 void evalTo(DesType& res) const 00618 { 00619 Index m = m_qr.rows(); 00620 Index n = m_qr.cols(); 00621 Index diagSize = (std::min)(m,n); 00622 res = m_other; 00623 if (m_transpose) 00624 { 00625 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 00626 //Compute res = Q' * other column by column 00627 for(Index j = 0; j < res.cols(); j++){ 00628 for (Index k = 0; k < diagSize; k++) 00629 { 00630 Scalar tau = Scalar(0); 00631 tau = m_qr.m_Q.col(k).dot(res.col(j)); 00632 if(tau==Scalar(0)) continue; 00633 tau = tau * m_qr.m_hcoeffs(k); 00634 res.col(j) -= tau * m_qr.m_Q.col(k); 00635 } 00636 } 00637 } 00638 else 00639 { 00640 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 00641 // Compute res = Q * other column by column 00642 for(Index j = 0; j < res.cols(); j++) 00643 { 00644 for (Index k = diagSize-1; k >=0; k--) 00645 { 00646 Scalar tau = Scalar(0); 00647 tau = m_qr.m_Q.col(k).dot(res.col(j)); 00648 if(tau==Scalar(0)) continue; 00649 tau = tau * m_qr.m_hcoeffs(k); 00650 res.col(j) -= tau * m_qr.m_Q.col(k); 00651 } 00652 } 00653 } 00654 } 00655 00656 const SparseQRType& m_qr; 00657 const Derived& m_other; 00658 bool m_transpose; 00659 }; 00660 00661 template<typename SparseQRType> 00662 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> > 00663 { 00664 typedef typename SparseQRType::Index Index; 00665 typedef typename SparseQRType::Scalar Scalar; 00666 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 00667 SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} 00668 template<typename Derived> 00669 SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) 00670 { 00671 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); 00672 } 00673 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const 00674 { 00675 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 00676 } 00677 inline Index rows() const { return m_qr.rows(); } 00678 inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } 00679 // To use for operations with the transpose of Q 00680 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const 00681 { 00682 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 00683 } 00684 template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const 00685 { 00686 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows()); 00687 } 00688 template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const 00689 { 00690 Dest idMat(m_qr.rows(), m_qr.rows()); 00691 idMat.setIdentity(); 00692 // Sort the sparse householder reflectors if needed 00693 const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q(); 00694 dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false); 00695 } 00696 00697 const SparseQRType& m_qr; 00698 }; 00699 00700 template<typename SparseQRType> 00701 struct SparseQRMatrixQTransposeReturnType 00702 { 00703 SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} 00704 template<typename Derived> 00705 SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other) 00706 { 00707 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true); 00708 } 00709 const SparseQRType& m_qr; 00710 }; 00711 00712 } // end namespace Eigen 00713 00714 #endif