ColPivHouseholderQR.h (25498B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr> 5 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com> 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_COLPIVOTINGHOUSEHOLDERQR_H 12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H 13 14 namespace Eigen { 15 16 namespace internal { 17 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> > 18 : traits<_MatrixType> 19 { 20 typedef MatrixXpr XprKind; 21 typedef SolverStorage StorageKind; 22 typedef int StorageIndex; 23 enum { Flags = 0 }; 24 }; 25 26 } // end namespace internal 27 28 /** \ingroup QR_Module 29 * 30 * \class ColPivHouseholderQR 31 * 32 * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting 33 * 34 * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition 35 * 36 * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R 37 * such that 38 * \f[ 39 * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R} 40 * \f] 41 * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 42 * upper triangular matrix. 43 * 44 * This decomposition performs column pivoting in order to be rank-revealing and improve 45 * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR. 46 * 47 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. 48 * 49 * \sa MatrixBase::colPivHouseholderQr() 50 */ 51 template<typename _MatrixType> class ColPivHouseholderQR 52 : public SolverBase<ColPivHouseholderQR<_MatrixType> > 53 { 54 public: 55 56 typedef _MatrixType MatrixType; 57 typedef SolverBase<ColPivHouseholderQR> Base; 58 friend class SolverBase<ColPivHouseholderQR>; 59 60 EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR) 61 enum { 62 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime, 63 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime 64 }; 65 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType; 66 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType; 67 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType; 68 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 69 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType; 70 typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType; 71 typedef typename MatrixType::PlainObject PlainObject; 72 73 private: 74 75 typedef typename PermutationType::StorageIndex PermIndexType; 76 77 public: 78 79 /** 80 * \brief Default Constructor. 81 * 82 * The default constructor is useful in cases in which the user intends to 83 * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&). 84 */ 85 ColPivHouseholderQR() 86 : m_qr(), 87 m_hCoeffs(), 88 m_colsPermutation(), 89 m_colsTranspositions(), 90 m_temp(), 91 m_colNormsUpdated(), 92 m_colNormsDirect(), 93 m_isInitialized(false), 94 m_usePrescribedThreshold(false) {} 95 96 /** \brief Default Constructor with memory preallocation 97 * 98 * Like the default constructor but with preallocation of the internal data 99 * according to the specified problem \a size. 100 * \sa ColPivHouseholderQR() 101 */ 102 ColPivHouseholderQR(Index rows, Index cols) 103 : m_qr(rows, cols), 104 m_hCoeffs((std::min)(rows,cols)), 105 m_colsPermutation(PermIndexType(cols)), 106 m_colsTranspositions(cols), 107 m_temp(cols), 108 m_colNormsUpdated(cols), 109 m_colNormsDirect(cols), 110 m_isInitialized(false), 111 m_usePrescribedThreshold(false) {} 112 113 /** \brief Constructs a QR factorization from a given matrix 114 * 115 * This constructor computes the QR factorization of the matrix \a matrix by calling 116 * the method compute(). It is a short cut for: 117 * 118 * \code 119 * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols()); 120 * qr.compute(matrix); 121 * \endcode 122 * 123 * \sa compute() 124 */ 125 template<typename InputType> 126 explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix) 127 : m_qr(matrix.rows(), matrix.cols()), 128 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 129 m_colsPermutation(PermIndexType(matrix.cols())), 130 m_colsTranspositions(matrix.cols()), 131 m_temp(matrix.cols()), 132 m_colNormsUpdated(matrix.cols()), 133 m_colNormsDirect(matrix.cols()), 134 m_isInitialized(false), 135 m_usePrescribedThreshold(false) 136 { 137 compute(matrix.derived()); 138 } 139 140 /** \brief Constructs a QR factorization from a given matrix 141 * 142 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. 143 * 144 * \sa ColPivHouseholderQR(const EigenBase&) 145 */ 146 template<typename InputType> 147 explicit ColPivHouseholderQR(EigenBase<InputType>& matrix) 148 : m_qr(matrix.derived()), 149 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())), 150 m_colsPermutation(PermIndexType(matrix.cols())), 151 m_colsTranspositions(matrix.cols()), 152 m_temp(matrix.cols()), 153 m_colNormsUpdated(matrix.cols()), 154 m_colNormsDirect(matrix.cols()), 155 m_isInitialized(false), 156 m_usePrescribedThreshold(false) 157 { 158 computeInPlace(); 159 } 160 161 #ifdef EIGEN_PARSED_BY_DOXYGEN 162 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which 163 * *this is the QR decomposition, if any exists. 164 * 165 * \param b the right-hand-side of the equation to solve. 166 * 167 * \returns a solution. 168 * 169 * \note_about_checking_solutions 170 * 171 * \note_about_arbitrary_choice_of_solution 172 * 173 * Example: \include ColPivHouseholderQR_solve.cpp 174 * Output: \verbinclude ColPivHouseholderQR_solve.out 175 */ 176 template<typename Rhs> 177 inline const Solve<ColPivHouseholderQR, Rhs> 178 solve(const MatrixBase<Rhs>& b) const; 179 #endif 180 181 HouseholderSequenceType householderQ() const; 182 HouseholderSequenceType matrixQ() const 183 { 184 return householderQ(); 185 } 186 187 /** \returns a reference to the matrix where the Householder QR decomposition is stored 188 */ 189 const MatrixType& matrixQR() const 190 { 191 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 192 return m_qr; 193 } 194 195 /** \returns a reference to the matrix where the result Householder QR is stored 196 * \warning The strict lower part of this matrix contains internal values. 197 * Only the upper triangular part should be referenced. To get it, use 198 * \code matrixR().template triangularView<Upper>() \endcode 199 * For rank-deficient matrices, use 200 * \code 201 * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 202 * \endcode 203 */ 204 const MatrixType& matrixR() const 205 { 206 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 207 return m_qr; 208 } 209 210 template<typename InputType> 211 ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix); 212 213 /** \returns a const reference to the column permutation matrix */ 214 const PermutationType& colsPermutation() const 215 { 216 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 217 return m_colsPermutation; 218 } 219 220 /** \returns the absolute value of the determinant of the matrix of which 221 * *this is the QR decomposition. It has only linear complexity 222 * (that is, O(n) where n is the dimension of the square matrix) 223 * as the QR decomposition has already been computed. 224 * 225 * \note This is only for square matrices. 226 * 227 * \warning a determinant can be very big or small, so for matrices 228 * of large enough dimension, there is a risk of overflow/underflow. 229 * One way to work around that is to use logAbsDeterminant() instead. 230 * 231 * \sa logAbsDeterminant(), MatrixBase::determinant() 232 */ 233 typename MatrixType::RealScalar absDeterminant() const; 234 235 /** \returns the natural log of the absolute value of the determinant of the matrix of which 236 * *this is the QR decomposition. It has only linear complexity 237 * (that is, O(n) where n is the dimension of the square matrix) 238 * as the QR decomposition has already been computed. 239 * 240 * \note This is only for square matrices. 241 * 242 * \note This method is useful to work around the risk of overflow/underflow that's inherent 243 * to determinant computation. 244 * 245 * \sa absDeterminant(), MatrixBase::determinant() 246 */ 247 typename MatrixType::RealScalar logAbsDeterminant() const; 248 249 /** \returns the rank of the matrix of which *this is the QR decomposition. 250 * 251 * \note This method has to determine which pivots should be considered nonzero. 252 * For that, it uses the threshold value that you can control by calling 253 * setThreshold(const RealScalar&). 254 */ 255 inline Index rank() const 256 { 257 using std::abs; 258 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 259 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold(); 260 Index result = 0; 261 for(Index i = 0; i < m_nonzero_pivots; ++i) 262 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold); 263 return result; 264 } 265 266 /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition. 267 * 268 * \note This method has to determine which pivots should be considered nonzero. 269 * For that, it uses the threshold value that you can control by calling 270 * setThreshold(const RealScalar&). 271 */ 272 inline Index dimensionOfKernel() const 273 { 274 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 275 return cols() - rank(); 276 } 277 278 /** \returns true if the matrix of which *this is the QR decomposition represents an injective 279 * linear map, i.e. has trivial kernel; false otherwise. 280 * 281 * \note This method has to determine which pivots should be considered nonzero. 282 * For that, it uses the threshold value that you can control by calling 283 * setThreshold(const RealScalar&). 284 */ 285 inline bool isInjective() const 286 { 287 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 288 return rank() == cols(); 289 } 290 291 /** \returns true if the matrix of which *this is the QR decomposition represents a surjective 292 * linear map; false otherwise. 293 * 294 * \note This method has to determine which pivots should be considered nonzero. 295 * For that, it uses the threshold value that you can control by calling 296 * setThreshold(const RealScalar&). 297 */ 298 inline bool isSurjective() const 299 { 300 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 301 return rank() == rows(); 302 } 303 304 /** \returns true if the matrix of which *this is the QR decomposition is invertible. 305 * 306 * \note This method has to determine which pivots should be considered nonzero. 307 * For that, it uses the threshold value that you can control by calling 308 * setThreshold(const RealScalar&). 309 */ 310 inline bool isInvertible() const 311 { 312 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 313 return isInjective() && isSurjective(); 314 } 315 316 /** \returns the inverse of the matrix of which *this is the QR decomposition. 317 * 318 * \note If this matrix is not invertible, the returned matrix has undefined coefficients. 319 * Use isInvertible() to first determine whether this matrix is invertible. 320 */ 321 inline const Inverse<ColPivHouseholderQR> inverse() const 322 { 323 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 324 return Inverse<ColPivHouseholderQR>(*this); 325 } 326 327 inline Index rows() const { return m_qr.rows(); } 328 inline Index cols() const { return m_qr.cols(); } 329 330 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q. 331 * 332 * For advanced uses only. 333 */ 334 const HCoeffsType& hCoeffs() const { return m_hCoeffs; } 335 336 /** Allows to prescribe a threshold to be used by certain methods, such as rank(), 337 * who need to determine when pivots are to be considered nonzero. This is not used for the 338 * QR decomposition itself. 339 * 340 * When it needs to get the threshold value, Eigen calls threshold(). By default, this 341 * uses a formula to automatically determine a reasonable threshold. 342 * Once you have called the present method setThreshold(const RealScalar&), 343 * your value is used instead. 344 * 345 * \param threshold The new value to use as the threshold. 346 * 347 * A pivot will be considered nonzero if its absolute value is strictly greater than 348 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ 349 * where maxpivot is the biggest pivot. 350 * 351 * If you want to come back to the default behavior, call setThreshold(Default_t) 352 */ 353 ColPivHouseholderQR& setThreshold(const RealScalar& threshold) 354 { 355 m_usePrescribedThreshold = true; 356 m_prescribedThreshold = threshold; 357 return *this; 358 } 359 360 /** Allows to come back to the default behavior, letting Eigen use its default formula for 361 * determining the threshold. 362 * 363 * You should pass the special object Eigen::Default as parameter here. 364 * \code qr.setThreshold(Eigen::Default); \endcode 365 * 366 * See the documentation of setThreshold(const RealScalar&). 367 */ 368 ColPivHouseholderQR& setThreshold(Default_t) 369 { 370 m_usePrescribedThreshold = false; 371 return *this; 372 } 373 374 /** Returns the threshold that will be used by certain methods such as rank(). 375 * 376 * See the documentation of setThreshold(const RealScalar&). 377 */ 378 RealScalar threshold() const 379 { 380 eigen_assert(m_isInitialized || m_usePrescribedThreshold); 381 return m_usePrescribedThreshold ? m_prescribedThreshold 382 // this formula comes from experimenting (see "LU precision tuning" thread on the list) 383 // and turns out to be identical to Higham's formula used already in LDLt. 384 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize()); 385 } 386 387 /** \returns the number of nonzero pivots in the QR decomposition. 388 * Here nonzero is meant in the exact sense, not in a fuzzy sense. 389 * So that notion isn't really intrinsically interesting, but it is 390 * still useful when implementing algorithms. 391 * 392 * \sa rank() 393 */ 394 inline Index nonzeroPivots() const 395 { 396 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 397 return m_nonzero_pivots; 398 } 399 400 /** \returns the absolute value of the biggest pivot, i.e. the biggest 401 * diagonal coefficient of R. 402 */ 403 RealScalar maxPivot() const { return m_maxpivot; } 404 405 /** \brief Reports whether the QR factorization was successful. 406 * 407 * \note This function always returns \c Success. It is provided for compatibility 408 * with other factorization routines. 409 * \returns \c Success 410 */ 411 ComputationInfo info() const 412 { 413 eigen_assert(m_isInitialized && "Decomposition is not initialized."); 414 return Success; 415 } 416 417 #ifndef EIGEN_PARSED_BY_DOXYGEN 418 template<typename RhsType, typename DstType> 419 void _solve_impl(const RhsType &rhs, DstType &dst) const; 420 421 template<bool Conjugate, typename RhsType, typename DstType> 422 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; 423 #endif 424 425 protected: 426 427 friend class CompleteOrthogonalDecomposition<MatrixType>; 428 429 static void check_template_parameters() 430 { 431 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 432 } 433 434 void computeInPlace(); 435 436 MatrixType m_qr; 437 HCoeffsType m_hCoeffs; 438 PermutationType m_colsPermutation; 439 IntRowVectorType m_colsTranspositions; 440 RowVectorType m_temp; 441 RealRowVectorType m_colNormsUpdated; 442 RealRowVectorType m_colNormsDirect; 443 bool m_isInitialized, m_usePrescribedThreshold; 444 RealScalar m_prescribedThreshold, m_maxpivot; 445 Index m_nonzero_pivots; 446 Index m_det_pq; 447 }; 448 449 template<typename MatrixType> 450 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const 451 { 452 using std::abs; 453 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 454 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 455 return abs(m_qr.diagonal().prod()); 456 } 457 458 template<typename MatrixType> 459 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const 460 { 461 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 462 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!"); 463 return m_qr.diagonal().cwiseAbs().array().log().sum(); 464 } 465 466 /** Performs the QR factorization of the given matrix \a matrix. The result of 467 * the factorization is stored into \c *this, and a reference to \c *this 468 * is returned. 469 * 470 * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&) 471 */ 472 template<typename MatrixType> 473 template<typename InputType> 474 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix) 475 { 476 m_qr = matrix.derived(); 477 computeInPlace(); 478 return *this; 479 } 480 481 template<typename MatrixType> 482 void ColPivHouseholderQR<MatrixType>::computeInPlace() 483 { 484 check_template_parameters(); 485 486 // the column permutation is stored as int indices, so just to be sure: 487 eigen_assert(m_qr.cols()<=NumTraits<int>::highest()); 488 489 using std::abs; 490 491 Index rows = m_qr.rows(); 492 Index cols = m_qr.cols(); 493 Index size = m_qr.diagonalSize(); 494 495 m_hCoeffs.resize(size); 496 497 m_temp.resize(cols); 498 499 m_colsTranspositions.resize(m_qr.cols()); 500 Index number_of_transpositions = 0; 501 502 m_colNormsUpdated.resize(cols); 503 m_colNormsDirect.resize(cols); 504 for (Index k = 0; k < cols; ++k) { 505 // colNormsDirect(k) caches the most recent directly computed norm of 506 // column k. 507 m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm(); 508 m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k); 509 } 510 511 RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows); 512 RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon()); 513 514 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case) 515 m_maxpivot = RealScalar(0); 516 517 for(Index k = 0; k < size; ++k) 518 { 519 // first, we look up in our table m_colNormsUpdated which column has the biggest norm 520 Index biggest_col_index; 521 RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index)); 522 biggest_col_index += k; 523 524 // Track the number of meaningful pivots but do not stop the decomposition to make 525 // sure that the initial matrix is properly reproduced. See bug 941. 526 if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k)) 527 m_nonzero_pivots = k; 528 529 // apply the transposition to the columns 530 m_colsTranspositions.coeffRef(k) = biggest_col_index; 531 if(k != biggest_col_index) { 532 m_qr.col(k).swap(m_qr.col(biggest_col_index)); 533 std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index)); 534 std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index)); 535 ++number_of_transpositions; 536 } 537 538 // generate the householder vector, store it below the diagonal 539 RealScalar beta; 540 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta); 541 542 // apply the householder transformation to the diagonal coefficient 543 m_qr.coeffRef(k,k) = beta; 544 545 // remember the maximum absolute value of diagonal coefficients 546 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta); 547 548 // apply the householder transformation 549 m_qr.bottomRightCorner(rows-k, cols-k-1) 550 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1)); 551 552 // update our table of norms of the columns 553 for (Index j = k + 1; j < cols; ++j) { 554 // The following implements the stable norm downgrade step discussed in 555 // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf 556 // and used in LAPACK routines xGEQPF and xGEQP3. 557 // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html 558 if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) { 559 RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j); 560 temp = (RealScalar(1) + temp) * (RealScalar(1) - temp); 561 temp = temp < RealScalar(0) ? RealScalar(0) : temp; 562 RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) / 563 m_colNormsDirect.coeffRef(j)); 564 if (temp2 <= norm_downdate_threshold) { 565 // The updated norm has become too inaccurate so re-compute the column 566 // norm directly. 567 m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm(); 568 m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j); 569 } else { 570 m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp); 571 } 572 } 573 } 574 } 575 576 m_colsPermutation.setIdentity(PermIndexType(cols)); 577 for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k) 578 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k))); 579 580 m_det_pq = (number_of_transpositions%2) ? -1 : 1; 581 m_isInitialized = true; 582 } 583 584 #ifndef EIGEN_PARSED_BY_DOXYGEN 585 template<typename _MatrixType> 586 template<typename RhsType, typename DstType> 587 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const 588 { 589 const Index nonzero_pivots = nonzeroPivots(); 590 591 if(nonzero_pivots == 0) 592 { 593 dst.setZero(); 594 return; 595 } 596 597 typename RhsType::PlainObject c(rhs); 598 599 c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() ); 600 601 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) 602 .template triangularView<Upper>() 603 .solveInPlace(c.topRows(nonzero_pivots)); 604 605 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i); 606 for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero(); 607 } 608 609 template<typename _MatrixType> 610 template<bool Conjugate, typename RhsType, typename DstType> 611 void ColPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const 612 { 613 const Index nonzero_pivots = nonzeroPivots(); 614 615 if(nonzero_pivots == 0) 616 { 617 dst.setZero(); 618 return; 619 } 620 621 typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs); 622 623 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots) 624 .template triangularView<Upper>() 625 .transpose().template conjugateIf<Conjugate>() 626 .solveInPlace(c.topRows(nonzero_pivots)); 627 628 dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots); 629 dst.bottomRows(rows()-nonzero_pivots).setZero(); 630 631 dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).template conjugateIf<!Conjugate>() ); 632 } 633 #endif 634 635 namespace internal { 636 637 template<typename DstXprType, typename MatrixType> 638 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense> 639 { 640 typedef ColPivHouseholderQR<MatrixType> QrType; 641 typedef Inverse<QrType> SrcXprType; 642 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &) 643 { 644 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols())); 645 } 646 }; 647 648 } // end namespace internal 649 650 /** \returns the matrix Q as a sequence of householder transformations. 651 * You can extract the meaningful part only by using: 652 * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/ 653 template<typename MatrixType> 654 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType> 655 ::householderQ() const 656 { 657 eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized."); 658 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate()); 659 } 660 661 /** \return the column-pivoting Householder QR decomposition of \c *this. 662 * 663 * \sa class ColPivHouseholderQR 664 */ 665 template<typename Derived> 666 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject> 667 MatrixBase<Derived>::colPivHouseholderQr() const 668 { 669 return ColPivHouseholderQR<PlainObject>(eval()); 670 } 671 672 } // end namespace Eigen 673 674 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H