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