![]() |
Eigen
3.3.3
|
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::StorageIndex StorageIndex; 00025 typedef typename ReturnType::StorageKind StorageKind; 00026 enum { 00027 RowsAtCompileTime = Dynamic, 00028 ColsAtCompileTime = Dynamic 00029 }; 00030 }; 00031 template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> > 00032 { 00033 typedef typename SparseQRType::MatrixType ReturnType; 00034 }; 00035 template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> > 00036 { 00037 typedef typename Derived::PlainObject ReturnType; 00038 }; 00039 } // End namespace internal 00040 00070 template<typename _MatrixType, typename _OrderingType> 00071 class SparseQR : public SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > 00072 { 00073 protected: 00074 typedef SparseSolverBase<SparseQR<_MatrixType,_OrderingType> > Base; 00075 using Base::m_isInitialized; 00076 public: 00077 using Base::_solve_impl; 00078 typedef _MatrixType MatrixType; 00079 typedef _OrderingType OrderingType; 00080 typedef typename MatrixType::Scalar Scalar; 00081 typedef typename MatrixType::RealScalar RealScalar; 00082 typedef typename MatrixType::StorageIndex StorageIndex; 00083 typedef SparseMatrix<Scalar,ColMajor,StorageIndex> QRMatrixType; 00084 typedef Matrix<StorageIndex, Dynamic, 1> IndexVector; 00085 typedef Matrix<Scalar, Dynamic, 1> ScalarVector; 00086 typedef PermutationMatrix<Dynamic, Dynamic, StorageIndex> PermutationType; 00087 00088 enum { 00089 ColsAtCompileTime = MatrixType::ColsAtCompileTime, 00090 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 00091 }; 00092 00093 public: 00094 SparseQR () : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) 00095 { } 00096 00103 explicit SparseQR(const MatrixType& mat) : m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false),m_isEtreeOk(false) 00104 { 00105 compute(mat); 00106 } 00107 00114 void compute(const MatrixType& mat) 00115 { 00116 analyzePattern(mat); 00117 factorize(mat); 00118 } 00119 void analyzePattern(const MatrixType& mat); 00120 void factorize(const MatrixType& mat); 00121 00124 inline Index rows() const { return m_pmat.rows(); } 00125 00128 inline Index cols() const { return m_pmat.cols();} 00129 00143 const QRMatrixType& matrixR() const { return m_R; } 00144 00149 Index rank() const 00150 { 00151 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00152 return m_nonzeropivots; 00153 } 00154 00173 SparseQRMatrixQReturnType<SparseQR> matrixQ() const 00174 { return SparseQRMatrixQReturnType<SparseQR>(*this); } 00175 00179 const PermutationType& colsPermutation() const 00180 { 00181 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00182 return m_outputPerm_c; 00183 } 00184 00188 std::string lastErrorMessage() const { return m_lastError; } 00189 00191 template<typename Rhs, typename Dest> 00192 bool _solve_impl(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const 00193 { 00194 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00195 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00196 00197 Index rank = this->rank(); 00198 00199 // Compute Q^T * b; 00200 typename Dest::PlainObject y, b; 00201 y = this->matrixQ().transpose() * B; 00202 b = y; 00203 00204 // Solve with the triangular matrix R 00205 y.resize((std::max<Index>)(cols(),y.rows()),y.cols()); 00206 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank)); 00207 y.bottomRows(y.rows()-rank).setZero(); 00208 00209 // Apply the column permutation 00210 if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); 00211 else dest = y.topRows(cols()); 00212 00213 m_info = Success; 00214 return true; 00215 } 00216 00222 void setPivotThreshold(const RealScalar& threshold) 00223 { 00224 m_useDefaultThreshold = false; 00225 m_threshold = threshold; 00226 } 00227 00232 template<typename Rhs> 00233 inline const Solve<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 00234 { 00235 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00236 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00237 return Solve<SparseQR, Rhs>(*this, B.derived()); 00238 } 00239 template<typename Rhs> 00240 inline const Solve<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 00241 { 00242 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 00243 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 00244 return Solve<SparseQR, Rhs>(*this, B.derived()); 00245 } 00246 00255 ComputationInfo info() const 00256 { 00257 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 00258 return m_info; 00259 } 00260 00261 00263 inline void _sort_matrix_Q() 00264 { 00265 if(this->m_isQSorted) return; 00266 // The matrix Q is sorted during the transposition 00267 SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q); 00268 this->m_Q = mQrm; 00269 this->m_isQSorted = true; 00270 } 00271 00272 00273 protected: 00274 bool m_analysisIsok; 00275 bool m_factorizationIsok; 00276 mutable ComputationInfo m_info; 00277 std::string m_lastError; 00278 QRMatrixType m_pmat; // Temporary matrix 00279 QRMatrixType m_R; // The triangular factor matrix 00280 QRMatrixType m_Q; // The orthogonal reflectors 00281 ScalarVector m_hcoeffs; // The Householder coefficients 00282 PermutationType m_perm_c; // Fill-reducing Column permutation 00283 PermutationType m_pivotperm; // The permutation for rank revealing 00284 PermutationType m_outputPerm_c; // The final column permutation 00285 RealScalar m_threshold; // Threshold to determine null Householder reflections 00286 bool m_useDefaultThreshold; // Use default threshold 00287 Index m_nonzeropivots; // Number of non zero pivots found 00288 IndexVector m_etree; // Column elimination tree 00289 IndexVector m_firstRowElt; // First element in each row 00290 bool m_isQSorted; // whether Q is sorted or not 00291 bool m_isEtreeOk; // whether the elimination tree match the initial input matrix 00292 00293 template <typename, typename > friend struct SparseQR_QProduct; 00294 00295 }; 00296 00306 template <typename MatrixType, typename OrderingType> 00307 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) 00308 { 00309 eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); 00310 // Copy to a column major matrix if the input is rowmajor 00311 typename internal::conditional<MatrixType::IsRowMajor,QRMatrixType,const MatrixType&>::type matCpy(mat); 00312 // Compute the column fill reducing ordering 00313 OrderingType ord; 00314 ord(matCpy, m_perm_c); 00315 Index n = mat.cols(); 00316 Index m = mat.rows(); 00317 Index diagSize = (std::min)(m,n); 00318 00319 if (!m_perm_c.size()) 00320 { 00321 m_perm_c.resize(n); 00322 m_perm_c.indices().setLinSpaced(n, 0,StorageIndex(n-1)); 00323 } 00324 00325 // Compute the column elimination tree of the permuted matrix 00326 m_outputPerm_c = m_perm_c.inverse(); 00327 internal::coletree(matCpy, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 00328 m_isEtreeOk = true; 00329 00330 m_R.resize(m, n); 00331 m_Q.resize(m, diagSize); 00332 00333 // Allocate space for nonzero elements : rough estimation 00334 m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree 00335 m_Q.reserve(2*mat.nonZeros()); 00336 m_hcoeffs.resize(diagSize); 00337 m_analysisIsok = true; 00338 } 00339 00347 template <typename MatrixType, typename OrderingType> 00348 void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) 00349 { 00350 using std::abs; 00351 00352 eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); 00353 StorageIndex m = StorageIndex(mat.rows()); 00354 StorageIndex n = StorageIndex(mat.cols()); 00355 StorageIndex diagSize = (std::min)(m,n); 00356 IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes 00357 IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q 00358 Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q 00359 ScalarVector tval(m); // The dense vector used to compute the current column 00360 RealScalar pivotThreshold = m_threshold; 00361 00362 m_R.setZero(); 00363 m_Q.setZero(); 00364 m_pmat = mat; 00365 if(!m_isEtreeOk) 00366 { 00367 m_outputPerm_c = m_perm_c.inverse(); 00368 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 00369 m_isEtreeOk = true; 00370 } 00371 00372 m_pmat.uncompress(); // To have the innerNonZeroPtr allocated 00373 00374 // Apply the fill-in reducing permutation lazily: 00375 { 00376 // If the input is row major, copy the original column indices, 00377 // otherwise directly use the input matrix 00378 // 00379 IndexVector originalOuterIndicesCpy; 00380 const StorageIndex *originalOuterIndices = mat.outerIndexPtr(); 00381 if(MatrixType::IsRowMajor) 00382 { 00383 originalOuterIndicesCpy = IndexVector::Map(m_pmat.outerIndexPtr(),n+1); 00384 originalOuterIndices = originalOuterIndicesCpy.data(); 00385 } 00386 00387 for (int i = 0; i < n; i++) 00388 { 00389 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; 00390 m_pmat.outerIndexPtr()[p] = originalOuterIndices[i]; 00391 m_pmat.innerNonZeroPtr()[p] = originalOuterIndices[i+1] - originalOuterIndices[i]; 00392 } 00393 } 00394 00395 /* Compute the default threshold as in MatLab, see: 00396 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing 00397 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 00398 */ 00399 if(m_useDefaultThreshold) 00400 { 00401 RealScalar max2Norm = 0.0; 00402 for (int j = 0; j < n; j++) max2Norm = numext::maxi(max2Norm, m_pmat.col(j).norm()); 00403 if(max2Norm==RealScalar(0)) 00404 max2Norm = RealScalar(1); 00405 pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); 00406 } 00407 00408 // Initialize the numerical permutation 00409 m_pivotperm.setIdentity(n); 00410 00411 StorageIndex nonzeroCol = 0; // Record the number of valid pivots 00412 m_Q.startVec(0); 00413 00414 // Left looking rank-revealing QR factorization: compute a column of R and Q at a time 00415 for (StorageIndex col = 0; col < n; ++col) 00416 { 00417 mark.setConstant(-1); 00418 m_R.startVec(col); 00419 mark(nonzeroCol) = col; 00420 Qidx(0) = nonzeroCol; 00421 nzcolR = 0; nzcolQ = 1; 00422 bool found_diag = nonzeroCol>=m; 00423 tval.setZero(); 00424 00425 // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., 00426 // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. 00427 // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, 00428 // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. 00429 for (typename QRMatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) 00430 { 00431 StorageIndex curIdx = nonzeroCol; 00432 if(itp) curIdx = StorageIndex(itp.row()); 00433 if(curIdx == nonzeroCol) found_diag = true; 00434 00435 // Get the nonzeros indexes of the current column of R 00436 StorageIndex st = m_firstRowElt(curIdx); // The traversal of the etree starts here 00437 if (st < 0 ) 00438 { 00439 m_lastError = "Empty row found during numerical factorization"; 00440 m_info = InvalidInput; 00441 return; 00442 } 00443 00444 // Traverse the etree 00445 Index bi = nzcolR; 00446 for (; mark(st) != col; st = m_etree(st)) 00447 { 00448 Ridx(nzcolR) = st; // Add this row to the list, 00449 mark(st) = col; // and mark this row as visited 00450 nzcolR++; 00451 } 00452 00453 // Reverse the list to get the topological ordering 00454 Index nt = nzcolR-bi; 00455 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1)); 00456 00457 // Copy the current (curIdx,pcol) value of the input matrix 00458 if(itp) tval(curIdx) = itp.value(); 00459 else tval(curIdx) = Scalar(0); 00460 00461 // Compute the pattern of Q(:,k) 00462 if(curIdx > nonzeroCol && mark(curIdx) != col ) 00463 { 00464 Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q, 00465 mark(curIdx) = col; // and mark it as visited 00466 nzcolQ++; 00467 } 00468 } 00469 00470 // Browse all the indexes of R(:,col) in reverse order 00471 for (Index i = nzcolR-1; i >= 0; i--) 00472 { 00473 Index curIdx = Ridx(i); 00474 00475 // Apply the curIdx-th householder vector to the current column (temporarily stored into tval) 00476 Scalar tdot(0); 00477 00478 // First compute q' * tval 00479 tdot = m_Q.col(curIdx).dot(tval); 00480 00481 tdot *= m_hcoeffs(curIdx); 00482 00483 // Then update tval = tval - q * tau 00484 // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse") 00485 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 00486 tval(itq.row()) -= itq.value() * tdot; 00487 00488 // Detect fill-in for the current column of Q 00489 if(m_etree(Ridx(i)) == nonzeroCol) 00490 { 00491 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 00492 { 00493 StorageIndex iQ = StorageIndex(itq.row()); 00494 if (mark(iQ) != col) 00495 { 00496 Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q, 00497 mark(iQ) = col; // and mark it as visited 00498 } 00499 } 00500 } 00501 } // End update current column 00502 00503 Scalar tau = RealScalar(0); 00504 RealScalar beta = 0; 00505 00506 if(nonzeroCol < diagSize) 00507 { 00508 // Compute the Householder reflection that eliminate the current column 00509 // FIXME this step should call the Householder module. 00510 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); 00511 00512 // First, the squared norm of Q((col+1):m, col) 00513 RealScalar sqrNorm = 0.; 00514 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); 00515 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) 00516 { 00517 beta = numext::real(c0); 00518 tval(Qidx(0)) = 1; 00519 } 00520 else 00521 { 00522 using std::sqrt; 00523 beta = sqrt(numext::abs2(c0) + sqrNorm); 00524 if(numext::real(c0) >= RealScalar(0)) 00525 beta = -beta; 00526 tval(Qidx(0)) = 1; 00527 for (Index itq = 1; itq < nzcolQ; ++itq) 00528 tval(Qidx(itq)) /= (c0 - beta); 00529 tau = numext::conj((beta-c0) / beta); 00530 00531 } 00532 } 00533 00534 // Insert values in R 00535 for (Index i = nzcolR-1; i >= 0; i--) 00536 { 00537 Index curIdx = Ridx(i); 00538 if(curIdx < nonzeroCol) 00539 { 00540 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); 00541 tval(curIdx) = Scalar(0.); 00542 } 00543 } 00544 00545 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold) 00546 { 00547 m_R.insertBackByOuterInner(col, nonzeroCol) = beta; 00548 // The householder coefficient 00549 m_hcoeffs(nonzeroCol) = tau; 00550 // Record the householder reflections 00551 for (Index itq = 0; itq < nzcolQ; ++itq) 00552 { 00553 Index iQ = Qidx(itq); 00554 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ); 00555 tval(iQ) = Scalar(0.); 00556 } 00557 nonzeroCol++; 00558 if(nonzeroCol<diagSize) 00559 m_Q.startVec(nonzeroCol); 00560 } 00561 else 00562 { 00563 // Zero pivot found: move implicitly this column to the end 00564 for (Index j = nonzeroCol; j < n-1; j++) 00565 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]); 00566 00567 // Recompute the column elimination tree 00568 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); 00569 m_isEtreeOk = false; 00570 } 00571 } 00572 00573 m_hcoeffs.tail(diagSize-nonzeroCol).setZero(); 00574 00575 // Finalize the column pointers of the sparse matrices R and Q 00576 m_Q.finalize(); 00577 m_Q.makeCompressed(); 00578 m_R.finalize(); 00579 m_R.makeCompressed(); 00580 m_isQSorted = false; 00581 00582 m_nonzeropivots = nonzeroCol; 00583 00584 if(nonzeroCol<n) 00585 { 00586 // Permute the triangular factor to put the 'dead' columns to the end 00587 QRMatrixType tempR(m_R); 00588 m_R = tempR * m_pivotperm; 00589 00590 // Update the column permutation 00591 m_outputPerm_c = m_outputPerm_c * m_pivotperm; 00592 } 00593 00594 m_isInitialized = true; 00595 m_factorizationIsok = true; 00596 m_info = Success; 00597 } 00598 00599 template <typename SparseQRType, typename Derived> 00600 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> > 00601 { 00602 typedef typename SparseQRType::QRMatrixType MatrixType; 00603 typedef typename SparseQRType::Scalar Scalar; 00604 // Get the references 00605 SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 00606 m_qr(qr),m_other(other),m_transpose(transpose) {} 00607 inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); } 00608 inline Index cols() const { return m_other.cols(); } 00609 00610 // Assign to a vector 00611 template<typename DesType> 00612 void evalTo(DesType& res) const 00613 { 00614 Index m = m_qr.rows(); 00615 Index n = m_qr.cols(); 00616 Index diagSize = (std::min)(m,n); 00617 res = m_other; 00618 if (m_transpose) 00619 { 00620 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 00621 //Compute res = Q' * other column by column 00622 for(Index j = 0; j < res.cols(); j++){ 00623 for (Index k = 0; k < diagSize; k++) 00624 { 00625 Scalar tau = Scalar(0); 00626 tau = m_qr.m_Q.col(k).dot(res.col(j)); 00627 if(tau==Scalar(0)) continue; 00628 tau = tau * m_qr.m_hcoeffs(k); 00629 res.col(j) -= tau * m_qr.m_Q.col(k); 00630 } 00631 } 00632 } 00633 else 00634 { 00635 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 00636 // Compute res = Q * other column by column 00637 for(Index j = 0; j < res.cols(); j++) 00638 { 00639 for (Index k = diagSize-1; k >=0; k--) 00640 { 00641 Scalar tau = Scalar(0); 00642 tau = m_qr.m_Q.col(k).dot(res.col(j)); 00643 if(tau==Scalar(0)) continue; 00644 tau = tau * m_qr.m_hcoeffs(k); 00645 res.col(j) -= tau * m_qr.m_Q.col(k); 00646 } 00647 } 00648 } 00649 } 00650 00651 const SparseQRType& m_qr; 00652 const Derived& m_other; 00653 bool m_transpose; 00654 }; 00655 00656 template<typename SparseQRType> 00657 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> > 00658 { 00659 typedef typename SparseQRType::Scalar Scalar; 00660 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 00661 enum { 00662 RowsAtCompileTime = Dynamic, 00663 ColsAtCompileTime = Dynamic 00664 }; 00665 explicit SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} 00666 template<typename Derived> 00667 SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) 00668 { 00669 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); 00670 } 00671 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const 00672 { 00673 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 00674 } 00675 inline Index rows() const { return m_qr.rows(); } 00676 inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } 00677 // To use for operations with the transpose of Q 00678 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const 00679 { 00680 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 00681 } 00682 const SparseQRType& m_qr; 00683 }; 00684 00685 template<typename SparseQRType> 00686 struct SparseQRMatrixQTransposeReturnType 00687 { 00688 explicit SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} 00689 template<typename Derived> 00690 SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other) 00691 { 00692 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true); 00693 } 00694 const SparseQRType& m_qr; 00695 }; 00696 00697 namespace internal { 00698 00699 template<typename SparseQRType> 00700 struct evaluator_traits<SparseQRMatrixQReturnType<SparseQRType> > 00701 { 00702 typedef typename SparseQRType::MatrixType MatrixType; 00703 typedef typename storage_kind_to_evaluator_kind<typename MatrixType::StorageKind>::Kind Kind; 00704 typedef SparseShape Shape; 00705 }; 00706 00707 template< typename DstXprType, typename SparseQRType> 00708 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Sparse> 00709 { 00710 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType; 00711 typedef typename DstXprType::Scalar Scalar; 00712 typedef typename DstXprType::StorageIndex StorageIndex; 00713 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/) 00714 { 00715 typename DstXprType::PlainObject idMat(src.m_qr.rows(), src.m_qr.rows()); 00716 idMat.setIdentity(); 00717 // Sort the sparse householder reflectors if needed 00718 const_cast<SparseQRType *>(&src.m_qr)->_sort_matrix_Q(); 00719 dst = SparseQR_QProduct<SparseQRType, DstXprType>(src.m_qr, idMat, false); 00720 } 00721 }; 00722 00723 template< typename DstXprType, typename SparseQRType> 00724 struct Assignment<DstXprType, SparseQRMatrixQReturnType<SparseQRType>, internal::assign_op<typename DstXprType::Scalar,typename DstXprType::Scalar>, Sparse2Dense> 00725 { 00726 typedef SparseQRMatrixQReturnType<SparseQRType> SrcXprType; 00727 typedef typename DstXprType::Scalar Scalar; 00728 typedef typename DstXprType::StorageIndex StorageIndex; 00729 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<Scalar,Scalar> &/*func*/) 00730 { 00731 dst = src.m_qr.matrixQ() * DstXprType::Identity(src.m_qr.rows(), src.m_qr.rows()); 00732 } 00733 }; 00734 00735 } // end namespace internal 00736 00737 } // end namespace Eigen 00738 00739 #endif