1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2012-2013 Desire Nuentsa <desire.nuentsa_wakam (at) inria.fr> 5 // Copyright (C) 2012-2014 Gael Guennebaud <gael.guennebaud (at) inria.fr> 6 // 7 // This Source Code Form is subject to the terms of the Mozilla 8 // Public License v. 2.0. If a copy of the MPL was not distributed 9 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 10 11 #ifndef EIGEN_SPARSE_QR_H 12 #define EIGEN_SPARSE_QR_H 13 14 namespace Eigen { 15 16 template<typename MatrixType, typename OrderingType> class SparseQR; 17 template<typename SparseQRType> struct SparseQRMatrixQReturnType; 18 template<typename SparseQRType> struct SparseQRMatrixQTransposeReturnType; 19 template<typename SparseQRType, typename Derived> struct SparseQR_QProduct; 20 namespace internal { 21 template <typename SparseQRType> struct traits<SparseQRMatrixQReturnType<SparseQRType> > 22 { 23 typedef typename SparseQRType::MatrixType ReturnType; 24 typedef typename ReturnType::Index Index; 25 typedef typename ReturnType::StorageKind StorageKind; 26 }; 27 template <typename SparseQRType> struct traits<SparseQRMatrixQTransposeReturnType<SparseQRType> > 28 { 29 typedef typename SparseQRType::MatrixType ReturnType; 30 }; 31 template <typename SparseQRType, typename Derived> struct traits<SparseQR_QProduct<SparseQRType, Derived> > 32 { 33 typedef typename Derived::PlainObject ReturnType; 34 }; 35 } // End namespace internal 36 37 /** 38 * \ingroup SparseQR_Module 39 * \class SparseQR 40 * \brief Sparse left-looking rank-revealing QR factorization 41 * 42 * This class implements a left-looking rank-revealing QR decomposition 43 * of sparse matrices. When a column has a norm less than a given tolerance 44 * it is implicitly permuted to the end. The QR factorization thus obtained is 45 * given by A*P = Q*R where R is upper triangular or trapezoidal. 46 * 47 * P is the column permutation which is the product of the fill-reducing and the 48 * rank-revealing permutations. Use colsPermutation() to get it. 49 * 50 * Q is the orthogonal matrix represented as products of Householder reflectors. 51 * Use matrixQ() to get an expression and matrixQ().transpose() to get the transpose. 52 * You can then apply it to a vector. 53 * 54 * R is the sparse triangular or trapezoidal matrix. The later occurs when A is rank-deficient. 55 * matrixR().topLeftCorner(rank(), rank()) always returns a triangular factor of full rank. 56 * 57 * \tparam _MatrixType The type of the sparse matrix A, must be a column-major SparseMatrix<> 58 * \tparam _OrderingType The fill-reducing ordering method. See the \link OrderingMethods_Module 59 * OrderingMethods \endlink module for the list of built-in and external ordering methods. 60 * 61 * \warning The input sparse matrix A must be in compressed mode (see SparseMatrix::makeCompressed()). 62 * 63 */ 64 template<typename _MatrixType, typename _OrderingType> 65 class SparseQR 66 { 67 public: 68 typedef _MatrixType MatrixType; 69 typedef _OrderingType OrderingType; 70 typedef typename MatrixType::Scalar Scalar; 71 typedef typename MatrixType::RealScalar RealScalar; 72 typedef typename MatrixType::Index Index; 73 typedef SparseMatrix<Scalar,ColMajor,Index> QRMatrixType; 74 typedef Matrix<Index, Dynamic, 1> IndexVector; 75 typedef Matrix<Scalar, Dynamic, 1> ScalarVector; 76 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType; 77 public: 78 SparseQR () : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) 79 { } 80 81 /** Construct a QR factorization of the matrix \a mat. 82 * 83 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 84 * 85 * \sa compute() 86 */ 87 SparseQR(const MatrixType& mat) : m_isInitialized(false), m_analysisIsok(false), m_lastError(""), m_useDefaultThreshold(true),m_isQSorted(false) 88 { 89 compute(mat); 90 } 91 92 /** Computes the QR factorization of the sparse matrix \a mat. 93 * 94 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 95 * 96 * \sa analyzePattern(), factorize() 97 */ 98 void compute(const MatrixType& mat) 99 { 100 analyzePattern(mat); 101 factorize(mat); 102 } 103 void analyzePattern(const MatrixType& mat); 104 void factorize(const MatrixType& mat); 105 106 /** \returns the number of rows of the represented matrix. 107 */ 108 inline Index rows() const { return m_pmat.rows(); } 109 110 /** \returns the number of columns of the represented matrix. 111 */ 112 inline Index cols() const { return m_pmat.cols();} 113 114 /** \returns a const reference to the \b sparse upper triangular matrix R of the QR factorization. 115 */ 116 const QRMatrixType& matrixR() const { return m_R; } 117 118 /** \returns the number of non linearly dependent columns as determined by the pivoting threshold. 119 * 120 * \sa setPivotThreshold() 121 */ 122 Index rank() const 123 { 124 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 125 return m_nonzeropivots; 126 } 127 128 /** \returns an expression of the matrix Q as products of sparse Householder reflectors. 129 * The common usage of this function is to apply it to a dense matrix or vector 130 * \code 131 * VectorXd B1, B2; 132 * // Initialize B1 133 * B2 = matrixQ() * B1; 134 * \endcode 135 * 136 * To get a plain SparseMatrix representation of Q: 137 * \code 138 * SparseMatrix<double> Q; 139 * Q = SparseQR<SparseMatrix<double> >(A).matrixQ(); 140 * \endcode 141 * Internally, this call simply performs a sparse product between the matrix Q 142 * and a sparse identity matrix. However, due to the fact that the sparse 143 * reflectors are stored unsorted, two transpositions are needed to sort 144 * them before performing the product. 145 */ 146 SparseQRMatrixQReturnType<SparseQR> matrixQ() const 147 { return SparseQRMatrixQReturnType<SparseQR>(*this); } 148 149 /** \returns a const reference to the column permutation P that was applied to A such that A*P = Q*R 150 * It is the combination of the fill-in reducing permutation and numerical column pivoting. 151 */ 152 const PermutationType& colsPermutation() const 153 { 154 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 155 return m_outputPerm_c; 156 } 157 158 /** \returns A string describing the type of error. 159 * This method is provided to ease debugging, not to handle errors. 160 */ 161 std::string lastErrorMessage() const { return m_lastError; } 162 163 /** \internal */ 164 template<typename Rhs, typename Dest> 165 bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &dest) const 166 { 167 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 168 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 169 170 Index rank = this->rank(); 171 172 // Compute Q^T * b; 173 typename Dest::PlainObject y, b; 174 y = this->matrixQ().transpose() * B; 175 b = y; 176 177 // Solve with the triangular matrix R 178 y.resize((std::max)(cols(),Index(y.rows())),y.cols()); 179 y.topRows(rank) = this->matrixR().topLeftCorner(rank, rank).template triangularView<Upper>().solve(b.topRows(rank)); 180 y.bottomRows(y.rows()-rank).setZero(); 181 182 // Apply the column permutation 183 if (m_perm_c.size()) dest = colsPermutation() * y.topRows(cols()); 184 else dest = y.topRows(cols()); 185 186 m_info = Success; 187 return true; 188 } 189 190 191 /** Sets the threshold that is used to determine linearly dependent columns during the factorization. 192 * 193 * In practice, if during the factorization the norm of the column that has to be eliminated is below 194 * this threshold, then the entire column is treated as zero, and it is moved at the end. 195 */ 196 void setPivotThreshold(const RealScalar& threshold) 197 { 198 m_useDefaultThreshold = false; 199 m_threshold = threshold; 200 } 201 202 /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A. 203 * 204 * \sa compute() 205 */ 206 template<typename Rhs> 207 inline const internal::solve_retval<SparseQR, Rhs> solve(const MatrixBase<Rhs>& B) const 208 { 209 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 210 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 211 return internal::solve_retval<SparseQR, Rhs>(*this, B.derived()); 212 } 213 template<typename Rhs> 214 inline const internal::sparse_solve_retval<SparseQR, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 215 { 216 eigen_assert(m_isInitialized && "The factorization should be called first, use compute()"); 217 eigen_assert(this->rows() == B.rows() && "SparseQR::solve() : invalid number of rows in the right hand side matrix"); 218 return internal::sparse_solve_retval<SparseQR, Rhs>(*this, B.derived()); 219 } 220 221 /** \brief Reports whether previous computation was successful. 222 * 223 * \returns \c Success if computation was successful, 224 * \c NumericalIssue if the QR factorization reports a numerical problem 225 * \c InvalidInput if the input matrix is invalid 226 * 227 * \sa iparm() 228 */ 229 ComputationInfo info() const 230 { 231 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 232 return m_info; 233 } 234 235 protected: 236 inline void sort_matrix_Q() 237 { 238 if(this->m_isQSorted) return; 239 // The matrix Q is sorted during the transposition 240 SparseMatrix<Scalar, RowMajor, Index> mQrm(this->m_Q); 241 this->m_Q = mQrm; 242 this->m_isQSorted = true; 243 } 244 245 246 protected: 247 bool m_isInitialized; 248 bool m_analysisIsok; 249 bool m_factorizationIsok; 250 mutable ComputationInfo m_info; 251 std::string m_lastError; 252 QRMatrixType m_pmat; // Temporary matrix 253 QRMatrixType m_R; // The triangular factor matrix 254 QRMatrixType m_Q; // The orthogonal reflectors 255 ScalarVector m_hcoeffs; // The Householder coefficients 256 PermutationType m_perm_c; // Fill-reducing Column permutation 257 PermutationType m_pivotperm; // The permutation for rank revealing 258 PermutationType m_outputPerm_c; // The final column permutation 259 RealScalar m_threshold; // Threshold to determine null Householder reflections 260 bool m_useDefaultThreshold; // Use default threshold 261 Index m_nonzeropivots; // Number of non zero pivots found 262 IndexVector m_etree; // Column elimination tree 263 IndexVector m_firstRowElt; // First element in each row 264 bool m_isQSorted; // whether Q is sorted or not 265 266 template <typename, typename > friend struct SparseQR_QProduct; 267 template <typename > friend struct SparseQRMatrixQReturnType; 268 269 }; 270 271 /** \brief Preprocessing step of a QR factorization 272 * 273 * \warning The matrix \a mat must be in compressed mode (see SparseMatrix::makeCompressed()). 274 * 275 * In this step, the fill-reducing permutation is computed and applied to the columns of A 276 * and the column elimination tree is computed as well. Only the sparsity pattern of \a mat is exploited. 277 * 278 * \note In this step it is assumed that there is no empty row in the matrix \a mat. 279 */ 280 template <typename MatrixType, typename OrderingType> 281 void SparseQR<MatrixType,OrderingType>::analyzePattern(const MatrixType& mat) 282 { 283 eigen_assert(mat.isCompressed() && "SparseQR requires a sparse matrix in compressed mode. Call .makeCompressed() before passing it to SparseQR"); 284 // Compute the column fill reducing ordering 285 OrderingType ord; 286 ord(mat, m_perm_c); 287 Index n = mat.cols(); 288 Index m = mat.rows(); 289 Index diagSize = (std::min)(m,n); 290 291 if (!m_perm_c.size()) 292 { 293 m_perm_c.resize(n); 294 m_perm_c.indices().setLinSpaced(n, 0,n-1); 295 } 296 297 // Compute the column elimination tree of the permuted matrix 298 m_outputPerm_c = m_perm_c.inverse(); 299 internal::coletree(mat, m_etree, m_firstRowElt, m_outputPerm_c.indices().data()); 300 301 m_R.resize(m, n); 302 m_Q.resize(m, diagSize); 303 304 // Allocate space for nonzero elements : rough estimation 305 m_R.reserve(2*mat.nonZeros()); //FIXME Get a more accurate estimation through symbolic factorization with the etree 306 m_Q.reserve(2*mat.nonZeros()); 307 m_hcoeffs.resize(diagSize); 308 m_analysisIsok = true; 309 } 310 311 /** \brief Performs the numerical QR factorization of the input matrix 312 * 313 * The function SparseQR::analyzePattern(const MatrixType&) must have been called beforehand with 314 * a matrix having the same sparsity pattern than \a mat. 315 * 316 * \param mat The sparse column-major matrix 317 */ 318 template <typename MatrixType, typename OrderingType> 319 void SparseQR<MatrixType,OrderingType>::factorize(const MatrixType& mat) 320 { 321 using std::abs; 322 using std::max; 323 324 eigen_assert(m_analysisIsok && "analyzePattern() should be called before this step"); 325 Index m = mat.rows(); 326 Index n = mat.cols(); 327 Index diagSize = (std::min)(m,n); 328 IndexVector mark((std::max)(m,n)); mark.setConstant(-1); // Record the visited nodes 329 IndexVector Ridx(n), Qidx(m); // Store temporarily the row indexes for the current column of R and Q 330 Index nzcolR, nzcolQ; // Number of nonzero for the current column of R and Q 331 ScalarVector tval(m); // The dense vector used to compute the current column 332 RealScalar pivotThreshold = m_threshold; 333 334 m_pmat = mat; 335 m_pmat.uncompress(); // To have the innerNonZeroPtr allocated 336 // Apply the fill-in reducing permutation lazily: 337 for (int i = 0; i < n; i++) 338 { 339 Index p = m_perm_c.size() ? m_perm_c.indices()(i) : i; 340 m_pmat.outerIndexPtr()[p] = mat.outerIndexPtr()[i]; 341 m_pmat.innerNonZeroPtr()[p] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 342 } 343 344 /* Compute the default threshold as in MatLab, see: 345 * Tim Davis, "Algorithm 915, SuiteSparseQR: Multifrontal Multithreaded Rank-Revealing 346 * Sparse QR Factorization, ACM Trans. on Math. Soft. 38(1), 2011, Page 8:3 347 */ 348 if(m_useDefaultThreshold) 349 { 350 RealScalar max2Norm = 0.0; 351 for (int j = 0; j < n; j++) max2Norm = (max)(max2Norm, m_pmat.col(j).norm()); 352 pivotThreshold = 20 * (m + n) * max2Norm * NumTraits<RealScalar>::epsilon(); 353 } 354 355 // Initialize the numerical permutation 356 m_pivotperm.setIdentity(n); 357 358 Index nonzeroCol = 0; // Record the number of valid pivots 359 m_Q.startVec(0); 360 361 // Left looking rank-revealing QR factorization: compute a column of R and Q at a time 362 for (Index col = 0; col < n; ++col) 363 { 364 mark.setConstant(-1); 365 m_R.startVec(col); 366 mark(nonzeroCol) = col; 367 Qidx(0) = nonzeroCol; 368 nzcolR = 0; nzcolQ = 1; 369 bool found_diag = nonzeroCol>=m; 370 tval.setZero(); 371 372 // Symbolic factorization: find the nonzero locations of the column k of the factors R and Q, i.e., 373 // all the nodes (with indexes lower than rank) reachable through the column elimination tree (etree) rooted at node k. 374 // Note: if the diagonal entry does not exist, then its contribution must be explicitly added, 375 // thus the trick with found_diag that permits to do one more iteration on the diagonal element if this one has not been found. 376 for (typename MatrixType::InnerIterator itp(m_pmat, col); itp || !found_diag; ++itp) 377 { 378 Index curIdx = nonzeroCol; 379 if(itp) curIdx = itp.row(); 380 if(curIdx == nonzeroCol) found_diag = true; 381 382 // Get the nonzeros indexes of the current column of R 383 Index st = m_firstRowElt(curIdx); // The traversal of the etree starts here 384 if (st < 0 ) 385 { 386 m_lastError = "Empty row found during numerical factorization"; 387 m_info = InvalidInput; 388 return; 389 } 390 391 // Traverse the etree 392 Index bi = nzcolR; 393 for (; mark(st) != col; st = m_etree(st)) 394 { 395 Ridx(nzcolR) = st; // Add this row to the list, 396 mark(st) = col; // and mark this row as visited 397 nzcolR++; 398 } 399 400 // Reverse the list to get the topological ordering 401 Index nt = nzcolR-bi; 402 for(Index i = 0; i < nt/2; i++) std::swap(Ridx(bi+i), Ridx(nzcolR-i-1)); 403 404 // Copy the current (curIdx,pcol) value of the input matrix 405 if(itp) tval(curIdx) = itp.value(); 406 else tval(curIdx) = Scalar(0); 407 408 // Compute the pattern of Q(:,k) 409 if(curIdx > nonzeroCol && mark(curIdx) != col ) 410 { 411 Qidx(nzcolQ) = curIdx; // Add this row to the pattern of Q, 412 mark(curIdx) = col; // and mark it as visited 413 nzcolQ++; 414 } 415 } 416 417 // Browse all the indexes of R(:,col) in reverse order 418 for (Index i = nzcolR-1; i >= 0; i--) 419 { 420 Index curIdx = Ridx(i); 421 422 // Apply the curIdx-th householder vector to the current column (temporarily stored into tval) 423 Scalar tdot(0); 424 425 // First compute q' * tval 426 tdot = m_Q.col(curIdx).dot(tval); 427 428 tdot *= m_hcoeffs(curIdx); 429 430 // Then update tval = tval - q * tau 431 // FIXME: tval -= tdot * m_Q.col(curIdx) should amount to the same (need to check/add support for efficient "dense ?= sparse") 432 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 433 tval(itq.row()) -= itq.value() * tdot; 434 435 // Detect fill-in for the current column of Q 436 if(m_etree(Ridx(i)) == nonzeroCol) 437 { 438 for (typename QRMatrixType::InnerIterator itq(m_Q, curIdx); itq; ++itq) 439 { 440 Index iQ = itq.row(); 441 if (mark(iQ) != col) 442 { 443 Qidx(nzcolQ++) = iQ; // Add this row to the pattern of Q, 444 mark(iQ) = col; // and mark it as visited 445 } 446 } 447 } 448 } // End update current column 449 450 Scalar tau; 451 RealScalar beta = 0; 452 453 if(nonzeroCol < diagSize) 454 { 455 // Compute the Householder reflection that eliminate the current column 456 // FIXME this step should call the Householder module. 457 Scalar c0 = nzcolQ ? tval(Qidx(0)) : Scalar(0); 458 459 // First, the squared norm of Q((col+1):m, col) 460 RealScalar sqrNorm = 0.; 461 for (Index itq = 1; itq < nzcolQ; ++itq) sqrNorm += numext::abs2(tval(Qidx(itq))); 462 if(sqrNorm == RealScalar(0) && numext::imag(c0) == RealScalar(0)) 463 { 464 tau = RealScalar(0); 465 beta = numext::real(c0); 466 tval(Qidx(0)) = 1; 467 } 468 else 469 { 470 using std::sqrt; 471 beta = sqrt(numext::abs2(c0) + sqrNorm); 472 if(numext::real(c0) >= RealScalar(0)) 473 beta = -beta; 474 tval(Qidx(0)) = 1; 475 for (Index itq = 1; itq < nzcolQ; ++itq) 476 tval(Qidx(itq)) /= (c0 - beta); 477 tau = numext::conj((beta-c0) / beta); 478 479 } 480 } 481 482 // Insert values in R 483 for (Index i = nzcolR-1; i >= 0; i--) 484 { 485 Index curIdx = Ridx(i); 486 if(curIdx < nonzeroCol) 487 { 488 m_R.insertBackByOuterInnerUnordered(col, curIdx) = tval(curIdx); 489 tval(curIdx) = Scalar(0.); 490 } 491 } 492 493 if(nonzeroCol < diagSize && abs(beta) >= pivotThreshold) 494 { 495 m_R.insertBackByOuterInner(col, nonzeroCol) = beta; 496 // The householder coefficient 497 m_hcoeffs(nonzeroCol) = tau; 498 // Record the householder reflections 499 for (Index itq = 0; itq < nzcolQ; ++itq) 500 { 501 Index iQ = Qidx(itq); 502 m_Q.insertBackByOuterInnerUnordered(nonzeroCol,iQ) = tval(iQ); 503 tval(iQ) = Scalar(0.); 504 } 505 nonzeroCol++; 506 if(nonzeroCol<diagSize) 507 m_Q.startVec(nonzeroCol); 508 } 509 else 510 { 511 // Zero pivot found: move implicitly this column to the end 512 for (Index j = nonzeroCol; j < n-1; j++) 513 std::swap(m_pivotperm.indices()(j), m_pivotperm.indices()[j+1]); 514 515 // Recompute the column elimination tree 516 internal::coletree(m_pmat, m_etree, m_firstRowElt, m_pivotperm.indices().data()); 517 } 518 } 519 520 m_hcoeffs.tail(diagSize-nonzeroCol).setZero(); 521 522 // Finalize the column pointers of the sparse matrices R and Q 523 m_Q.finalize(); 524 m_Q.makeCompressed(); 525 m_R.finalize(); 526 m_R.makeCompressed(); 527 m_isQSorted = false; 528 529 m_nonzeropivots = nonzeroCol; 530 531 if(nonzeroCol<n) 532 { 533 // Permute the triangular factor to put the 'dead' columns to the end 534 MatrixType tempR(m_R); 535 m_R = tempR * m_pivotperm; 536 537 // Update the column permutation 538 m_outputPerm_c = m_outputPerm_c * m_pivotperm; 539 } 540 541 m_isInitialized = true; 542 m_factorizationIsok = true; 543 m_info = Success; 544 } 545 546 namespace internal { 547 548 template<typename _MatrixType, typename OrderingType, typename Rhs> 549 struct solve_retval<SparseQR<_MatrixType,OrderingType>, Rhs> 550 : solve_retval_base<SparseQR<_MatrixType,OrderingType>, Rhs> 551 { 552 typedef SparseQR<_MatrixType,OrderingType> Dec; 553 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 554 555 template<typename Dest> void evalTo(Dest& dst) const 556 { 557 dec()._solve(rhs(),dst); 558 } 559 }; 560 template<typename _MatrixType, typename OrderingType, typename Rhs> 561 struct sparse_solve_retval<SparseQR<_MatrixType, OrderingType>, Rhs> 562 : sparse_solve_retval_base<SparseQR<_MatrixType, OrderingType>, Rhs> 563 { 564 typedef SparseQR<_MatrixType, OrderingType> Dec; 565 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec, Rhs) 566 567 template<typename Dest> void evalTo(Dest& dst) const 568 { 569 this->defaultEvalTo(dst); 570 } 571 }; 572 } // end namespace internal 573 574 template <typename SparseQRType, typename Derived> 575 struct SparseQR_QProduct : ReturnByValue<SparseQR_QProduct<SparseQRType, Derived> > 576 { 577 typedef typename SparseQRType::QRMatrixType MatrixType; 578 typedef typename SparseQRType::Scalar Scalar; 579 typedef typename SparseQRType::Index Index; 580 // Get the references 581 SparseQR_QProduct(const SparseQRType& qr, const Derived& other, bool transpose) : 582 m_qr(qr),m_other(other),m_transpose(transpose) {} 583 inline Index rows() const { return m_transpose ? m_qr.rows() : m_qr.cols(); } 584 inline Index cols() const { return m_other.cols(); } 585 586 // Assign to a vector 587 template<typename DesType> 588 void evalTo(DesType& res) const 589 { 590 Index m = m_qr.rows(); 591 Index n = m_qr.cols(); 592 Index diagSize = (std::min)(m,n); 593 res = m_other; 594 if (m_transpose) 595 { 596 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 597 //Compute res = Q' * other column by column 598 for(Index j = 0; j < res.cols(); j++){ 599 for (Index k = 0; k < diagSize; k++) 600 { 601 Scalar tau = Scalar(0); 602 tau = m_qr.m_Q.col(k).dot(res.col(j)); 603 if(tau==Scalar(0)) continue; 604 tau = tau * m_qr.m_hcoeffs(k); 605 res.col(j) -= tau * m_qr.m_Q.col(k); 606 } 607 } 608 } 609 else 610 { 611 eigen_assert(m_qr.m_Q.rows() == m_other.rows() && "Non conforming object sizes"); 612 // Compute res = Q * other column by column 613 for(Index j = 0; j < res.cols(); j++) 614 { 615 for (Index k = diagSize-1; k >=0; k--) 616 { 617 Scalar tau = Scalar(0); 618 tau = m_qr.m_Q.col(k).dot(res.col(j)); 619 if(tau==Scalar(0)) continue; 620 tau = tau * m_qr.m_hcoeffs(k); 621 res.col(j) -= tau * m_qr.m_Q.col(k); 622 } 623 } 624 } 625 } 626 627 const SparseQRType& m_qr; 628 const Derived& m_other; 629 bool m_transpose; 630 }; 631 632 template<typename SparseQRType> 633 struct SparseQRMatrixQReturnType : public EigenBase<SparseQRMatrixQReturnType<SparseQRType> > 634 { 635 typedef typename SparseQRType::Index Index; 636 typedef typename SparseQRType::Scalar Scalar; 637 typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix; 638 SparseQRMatrixQReturnType(const SparseQRType& qr) : m_qr(qr) {} 639 template<typename Derived> 640 SparseQR_QProduct<SparseQRType, Derived> operator*(const MatrixBase<Derived>& other) 641 { 642 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(),false); 643 } 644 SparseQRMatrixQTransposeReturnType<SparseQRType> adjoint() const 645 { 646 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 647 } 648 inline Index rows() const { return m_qr.rows(); } 649 inline Index cols() const { return (std::min)(m_qr.rows(),m_qr.cols()); } 650 // To use for operations with the transpose of Q 651 SparseQRMatrixQTransposeReturnType<SparseQRType> transpose() const 652 { 653 return SparseQRMatrixQTransposeReturnType<SparseQRType>(m_qr); 654 } 655 template<typename Dest> void evalTo(MatrixBase<Dest>& dest) const 656 { 657 dest.derived() = m_qr.matrixQ() * Dest::Identity(m_qr.rows(), m_qr.rows()); 658 } 659 template<typename Dest> void evalTo(SparseMatrixBase<Dest>& dest) const 660 { 661 Dest idMat(m_qr.rows(), m_qr.rows()); 662 idMat.setIdentity(); 663 // Sort the sparse householder reflectors if needed 664 const_cast<SparseQRType *>(&m_qr)->sort_matrix_Q(); 665 dest.derived() = SparseQR_QProduct<SparseQRType, Dest>(m_qr, idMat, false); 666 } 667 668 const SparseQRType& m_qr; 669 }; 670 671 template<typename SparseQRType> 672 struct SparseQRMatrixQTransposeReturnType 673 { 674 SparseQRMatrixQTransposeReturnType(const SparseQRType& qr) : m_qr(qr) {} 675 template<typename Derived> 676 SparseQR_QProduct<SparseQRType,Derived> operator*(const MatrixBase<Derived>& other) 677 { 678 return SparseQR_QProduct<SparseQRType,Derived>(m_qr,other.derived(), true); 679 } 680 const SparseQRType& m_qr; 681 }; 682 683 } // end namespace Eigen 684 685 #endif 686