CompleteOrthogonalDecomposition.h (23429B)
1 // This file is part of Eigen, a lightweight C++ template library 2 // for linear algebra. 3 // 4 // Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com> 5 // 6 // This Source Code Form is subject to the terms of the Mozilla 7 // Public License v. 2.0. If a copy of the MPL was not distributed 8 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 9 10 #ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H 11 #define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H 12 13 namespace Eigen { 14 15 namespace internal { 16 template <typename _MatrixType> 17 struct traits<CompleteOrthogonalDecomposition<_MatrixType> > 18 : traits<_MatrixType> { 19 typedef MatrixXpr XprKind; 20 typedef SolverStorage StorageKind; 21 typedef int StorageIndex; 22 enum { Flags = 0 }; 23 }; 24 25 } // end namespace internal 26 27 /** \ingroup QR_Module 28 * 29 * \class CompleteOrthogonalDecomposition 30 * 31 * \brief Complete orthogonal decomposition (COD) of a matrix. 32 * 33 * \param MatrixType the type of the matrix of which we are computing the COD. 34 * 35 * This class performs a rank-revealing complete orthogonal decomposition of a 36 * matrix \b A into matrices \b P, \b Q, \b T, and \b Z such that 37 * \f[ 38 * \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, 39 * \begin{bmatrix} \mathbf{T} & \mathbf{0} \\ 40 * \mathbf{0} & \mathbf{0} \end{bmatrix} \, \mathbf{Z} 41 * \f] 42 * by using Householder transformations. Here, \b P is a permutation matrix, 43 * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of 44 * size rank-by-rank. \b A may be rank deficient. 45 * 46 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism. 47 * 48 * \sa MatrixBase::completeOrthogonalDecomposition() 49 */ 50 template <typename _MatrixType> class CompleteOrthogonalDecomposition 51 : public SolverBase<CompleteOrthogonalDecomposition<_MatrixType> > 52 { 53 public: 54 typedef _MatrixType MatrixType; 55 typedef SolverBase<CompleteOrthogonalDecomposition> Base; 56 57 template<typename Derived> 58 friend struct internal::solve_assertion; 59 60 EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition) 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> 67 PermutationType; 68 typedef typename internal::plain_row_type<MatrixType, Index>::type 69 IntRowVectorType; 70 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType; 71 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type 72 RealRowVectorType; 73 typedef HouseholderSequence< 74 MatrixType, typename internal::remove_all< 75 typename HCoeffsType::ConjugateReturnType>::type> 76 HouseholderSequenceType; 77 typedef typename MatrixType::PlainObject PlainObject; 78 79 private: 80 typedef typename PermutationType::Index PermIndexType; 81 82 public: 83 /** 84 * \brief Default Constructor. 85 * 86 * The default constructor is useful in cases in which the user intends to 87 * perform decompositions via 88 * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&). 89 */ 90 CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {} 91 92 /** \brief Default Constructor with memory preallocation 93 * 94 * Like the default constructor but with preallocation of the internal data 95 * according to the specified problem \a size. 96 * \sa CompleteOrthogonalDecomposition() 97 */ 98 CompleteOrthogonalDecomposition(Index rows, Index cols) 99 : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {} 100 101 /** \brief Constructs a complete orthogonal decomposition from a given 102 * matrix. 103 * 104 * This constructor computes the complete orthogonal decomposition of the 105 * matrix \a matrix by calling the method compute(). The default 106 * threshold for rank determination will be used. It is a short cut for: 107 * 108 * \code 109 * CompleteOrthogonalDecomposition<MatrixType> cod(matrix.rows(), 110 * matrix.cols()); 111 * cod.setThreshold(Default); 112 * cod.compute(matrix); 113 * \endcode 114 * 115 * \sa compute() 116 */ 117 template <typename InputType> 118 explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix) 119 : m_cpqr(matrix.rows(), matrix.cols()), 120 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), 121 m_temp(matrix.cols()) 122 { 123 compute(matrix.derived()); 124 } 125 126 /** \brief Constructs a complete orthogonal decomposition from a given matrix 127 * 128 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref. 129 * 130 * \sa CompleteOrthogonalDecomposition(const EigenBase&) 131 */ 132 template<typename InputType> 133 explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix) 134 : m_cpqr(matrix.derived()), 135 m_zCoeffs((std::min)(matrix.rows(), matrix.cols())), 136 m_temp(matrix.cols()) 137 { 138 computeInPlace(); 139 } 140 141 #ifdef EIGEN_PARSED_BY_DOXYGEN 142 /** This method computes the minimum-norm solution X to a least squares 143 * problem \f[\mathrm{minimize} \|A X - B\|, \f] where \b A is the matrix of 144 * which \c *this is the complete orthogonal decomposition. 145 * 146 * \param b the right-hand sides of the problem to solve. 147 * 148 * \returns a solution. 149 * 150 */ 151 template <typename Rhs> 152 inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve( 153 const MatrixBase<Rhs>& b) const; 154 #endif 155 156 HouseholderSequenceType householderQ(void) const; 157 HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); } 158 159 /** \returns the matrix \b Z. 160 */ 161 MatrixType matrixZ() const { 162 MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols()); 163 applyZOnTheLeftInPlace<false>(Z); 164 return Z; 165 } 166 167 /** \returns a reference to the matrix where the complete orthogonal 168 * decomposition is stored 169 */ 170 const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); } 171 172 /** \returns a reference to the matrix where the complete orthogonal 173 * decomposition is stored. 174 * \warning The strict lower part and \code cols() - rank() \endcode right 175 * columns of this matrix contains internal values. 176 * Only the upper triangular part should be referenced. To get it, use 177 * \code matrixT().template triangularView<Upper>() \endcode 178 * For rank-deficient matrices, use 179 * \code 180 * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 181 * \endcode 182 */ 183 const MatrixType& matrixT() const { return m_cpqr.matrixQR(); } 184 185 template <typename InputType> 186 CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) { 187 // Compute the column pivoted QR factorization A P = Q R. 188 m_cpqr.compute(matrix); 189 computeInPlace(); 190 return *this; 191 } 192 193 /** \returns a const reference to the column permutation matrix */ 194 const PermutationType& colsPermutation() const { 195 return m_cpqr.colsPermutation(); 196 } 197 198 /** \returns the absolute value of the determinant of the matrix of which 199 * *this is the complete orthogonal decomposition. It has only linear 200 * complexity (that is, O(n) where n is the dimension of the square matrix) 201 * as the complete orthogonal decomposition has already been computed. 202 * 203 * \note This is only for square matrices. 204 * 205 * \warning a determinant can be very big or small, so for matrices 206 * of large enough dimension, there is a risk of overflow/underflow. 207 * One way to work around that is to use logAbsDeterminant() instead. 208 * 209 * \sa logAbsDeterminant(), MatrixBase::determinant() 210 */ 211 typename MatrixType::RealScalar absDeterminant() const; 212 213 /** \returns the natural log of the absolute value of the determinant of the 214 * matrix of which *this is the complete orthogonal decomposition. It has 215 * only linear complexity (that is, O(n) where n is the dimension of the 216 * square matrix) as the complete orthogonal decomposition has already been 217 * computed. 218 * 219 * \note This is only for square matrices. 220 * 221 * \note This method is useful to work around the risk of overflow/underflow 222 * that's inherent to determinant computation. 223 * 224 * \sa absDeterminant(), MatrixBase::determinant() 225 */ 226 typename MatrixType::RealScalar logAbsDeterminant() const; 227 228 /** \returns the rank of the matrix of which *this is the complete orthogonal 229 * decomposition. 230 * 231 * \note This method has to determine which pivots should be considered 232 * nonzero. For that, it uses the threshold value that you can control by 233 * calling setThreshold(const RealScalar&). 234 */ 235 inline Index rank() const { return m_cpqr.rank(); } 236 237 /** \returns the dimension of the kernel of the matrix of which *this is the 238 * complete orthogonal decomposition. 239 * 240 * \note This method has to determine which pivots should be considered 241 * nonzero. For that, it uses the threshold value that you can control by 242 * calling setThreshold(const RealScalar&). 243 */ 244 inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); } 245 246 /** \returns true if the matrix of which *this is the decomposition represents 247 * an injective linear map, i.e. has trivial kernel; false otherwise. 248 * 249 * \note This method has to determine which pivots should be considered 250 * nonzero. For that, it uses the threshold value that you can control by 251 * calling setThreshold(const RealScalar&). 252 */ 253 inline bool isInjective() const { return m_cpqr.isInjective(); } 254 255 /** \returns true if the matrix of which *this is the decomposition represents 256 * a surjective linear map; false otherwise. 257 * 258 * \note This method has to determine which pivots should be considered 259 * nonzero. For that, it uses the threshold value that you can control by 260 * calling setThreshold(const RealScalar&). 261 */ 262 inline bool isSurjective() const { return m_cpqr.isSurjective(); } 263 264 /** \returns true if the matrix of which *this is the complete orthogonal 265 * decomposition is invertible. 266 * 267 * \note This method has to determine which pivots should be considered 268 * nonzero. For that, it uses the threshold value that you can control by 269 * calling setThreshold(const RealScalar&). 270 */ 271 inline bool isInvertible() const { return m_cpqr.isInvertible(); } 272 273 /** \returns the pseudo-inverse of the matrix of which *this is the complete 274 * orthogonal decomposition. 275 * \warning: Do not compute \c this->pseudoInverse()*rhs to solve a linear systems. 276 * It is more efficient and numerically stable to call \c this->solve(rhs). 277 */ 278 inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const 279 { 280 eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized."); 281 return Inverse<CompleteOrthogonalDecomposition>(*this); 282 } 283 284 inline Index rows() const { return m_cpqr.rows(); } 285 inline Index cols() const { return m_cpqr.cols(); } 286 287 /** \returns a const reference to the vector of Householder coefficients used 288 * to represent the factor \c Q. 289 * 290 * For advanced uses only. 291 */ 292 inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); } 293 294 /** \returns a const reference to the vector of Householder coefficients 295 * used to represent the factor \c Z. 296 * 297 * For advanced uses only. 298 */ 299 const HCoeffsType& zCoeffs() const { return m_zCoeffs; } 300 301 /** Allows to prescribe a threshold to be used by certain methods, such as 302 * rank(), who need to determine when pivots are to be considered nonzero. 303 * Most be called before calling compute(). 304 * 305 * When it needs to get the threshold value, Eigen calls threshold(). By 306 * default, this uses a formula to automatically determine a reasonable 307 * threshold. Once you have called the present method 308 * setThreshold(const RealScalar&), your value is used instead. 309 * 310 * \param threshold The new value to use as the threshold. 311 * 312 * A pivot will be considered nonzero if its absolute value is strictly 313 * greater than 314 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$ 315 * where maxpivot is the biggest pivot. 316 * 317 * If you want to come back to the default behavior, call 318 * setThreshold(Default_t) 319 */ 320 CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) { 321 m_cpqr.setThreshold(threshold); 322 return *this; 323 } 324 325 /** Allows to come back to the default behavior, letting Eigen use its default 326 * formula for determining the threshold. 327 * 328 * You should pass the special object Eigen::Default as parameter here. 329 * \code qr.setThreshold(Eigen::Default); \endcode 330 * 331 * See the documentation of setThreshold(const RealScalar&). 332 */ 333 CompleteOrthogonalDecomposition& setThreshold(Default_t) { 334 m_cpqr.setThreshold(Default); 335 return *this; 336 } 337 338 /** Returns the threshold that will be used by certain methods such as rank(). 339 * 340 * See the documentation of setThreshold(const RealScalar&). 341 */ 342 RealScalar threshold() const { return m_cpqr.threshold(); } 343 344 /** \returns the number of nonzero pivots in the complete orthogonal 345 * decomposition. Here nonzero is meant in the exact sense, not in a 346 * fuzzy sense. So that notion isn't really intrinsically interesting, 347 * but it is still useful when implementing algorithms. 348 * 349 * \sa rank() 350 */ 351 inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); } 352 353 /** \returns the absolute value of the biggest pivot, i.e. the biggest 354 * diagonal coefficient of R. 355 */ 356 inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); } 357 358 /** \brief Reports whether the complete orthogonal decomposition was 359 * successful. 360 * 361 * \note This function always returns \c Success. It is provided for 362 * compatibility 363 * with other factorization routines. 364 * \returns \c Success 365 */ 366 ComputationInfo info() const { 367 eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized."); 368 return Success; 369 } 370 371 #ifndef EIGEN_PARSED_BY_DOXYGEN 372 template <typename RhsType, typename DstType> 373 void _solve_impl(const RhsType& rhs, DstType& dst) const; 374 375 template<bool Conjugate, typename RhsType, typename DstType> 376 void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const; 377 #endif 378 379 protected: 380 static void check_template_parameters() { 381 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar); 382 } 383 384 template<bool Transpose_, typename Rhs> 385 void _check_solve_assertion(const Rhs& b) const { 386 EIGEN_ONLY_USED_FOR_DEBUG(b); 387 eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized."); 388 eigen_assert((Transpose_?derived().cols():derived().rows())==b.rows() && "CompleteOrthogonalDecomposition::solve(): invalid number of rows of the right hand side matrix b"); 389 } 390 391 void computeInPlace(); 392 393 /** Overwrites \b rhs with \f$ \mathbf{Z} * \mathbf{rhs} \f$ or 394 * \f$ \mathbf{\overline Z} * \mathbf{rhs} \f$ if \c Conjugate 395 * is set to \c true. 396 */ 397 template <bool Conjugate, typename Rhs> 398 void applyZOnTheLeftInPlace(Rhs& rhs) const; 399 400 /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$. 401 */ 402 template <typename Rhs> 403 void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const; 404 405 ColPivHouseholderQR<MatrixType> m_cpqr; 406 HCoeffsType m_zCoeffs; 407 RowVectorType m_temp; 408 }; 409 410 template <typename MatrixType> 411 typename MatrixType::RealScalar 412 CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const { 413 return m_cpqr.absDeterminant(); 414 } 415 416 template <typename MatrixType> 417 typename MatrixType::RealScalar 418 CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const { 419 return m_cpqr.logAbsDeterminant(); 420 } 421 422 /** Performs the complete orthogonal decomposition of the given matrix \a 423 * matrix. The result of the factorization is stored into \c *this, and a 424 * reference to \c *this is returned. 425 * 426 * \sa class CompleteOrthogonalDecomposition, 427 * CompleteOrthogonalDecomposition(const MatrixType&) 428 */ 429 template <typename MatrixType> 430 void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace() 431 { 432 check_template_parameters(); 433 434 // the column permutation is stored as int indices, so just to be sure: 435 eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest()); 436 437 const Index rank = m_cpqr.rank(); 438 const Index cols = m_cpqr.cols(); 439 const Index rows = m_cpqr.rows(); 440 m_zCoeffs.resize((std::min)(rows, cols)); 441 m_temp.resize(cols); 442 443 if (rank < cols) { 444 // We have reduced the (permuted) matrix to the form 445 // [R11 R12] 446 // [ 0 R22] 447 // where R11 is r-by-r (r = rank) upper triangular, R12 is 448 // r-by-(n-r), and R22 is empty or the norm of R22 is negligible. 449 // We now compute the complete orthogonal decomposition by applying 450 // Householder transformations from the right to the upper trapezoidal 451 // matrix X = [R11 R12] to zero out R12 and obtain the factorization 452 // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and 453 // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix. 454 // We store the data representing Z in R12 and m_zCoeffs. 455 for (Index k = rank - 1; k >= 0; --k) { 456 if (k != rank - 1) { 457 // Given the API for Householder reflectors, it is more convenient if 458 // we swap the leading parts of columns k and r-1 (zero-based) to form 459 // the matrix X_k = [X(0:k, k), X(0:k, r:n)] 460 m_cpqr.m_qr.col(k).head(k + 1).swap( 461 m_cpqr.m_qr.col(rank - 1).head(k + 1)); 462 } 463 // Construct Householder reflector Z(k) to zero out the last row of X_k, 464 // i.e. choose Z(k) such that 465 // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0]. 466 RealScalar beta; 467 m_cpqr.m_qr.row(k) 468 .tail(cols - rank + 1) 469 .makeHouseholderInPlace(m_zCoeffs(k), beta); 470 m_cpqr.m_qr(k, rank - 1) = beta; 471 if (k > 0) { 472 // Apply Z(k) to the first k rows of X_k 473 m_cpqr.m_qr.topRightCorner(k, cols - rank + 1) 474 .applyHouseholderOnTheRight( 475 m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k), 476 &m_temp(0)); 477 } 478 if (k != rank - 1) { 479 // Swap X(0:k,k) back to its proper location. 480 m_cpqr.m_qr.col(k).head(k + 1).swap( 481 m_cpqr.m_qr.col(rank - 1).head(k + 1)); 482 } 483 } 484 } 485 } 486 487 template <typename MatrixType> 488 template <bool Conjugate, typename Rhs> 489 void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace( 490 Rhs& rhs) const { 491 const Index cols = this->cols(); 492 const Index nrhs = rhs.cols(); 493 const Index rank = this->rank(); 494 Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs)); 495 for (Index k = rank-1; k >= 0; --k) { 496 if (k != rank - 1) { 497 rhs.row(k).swap(rhs.row(rank - 1)); 498 } 499 rhs.middleRows(rank - 1, cols - rank + 1) 500 .applyHouseholderOnTheLeft( 501 matrixQTZ().row(k).tail(cols - rank).transpose().template conjugateIf<!Conjugate>(), zCoeffs().template conjugateIf<Conjugate>()(k), 502 &temp(0)); 503 if (k != rank - 1) { 504 rhs.row(k).swap(rhs.row(rank - 1)); 505 } 506 } 507 } 508 509 template <typename MatrixType> 510 template <typename Rhs> 511 void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace( 512 Rhs& rhs) const { 513 const Index cols = this->cols(); 514 const Index nrhs = rhs.cols(); 515 const Index rank = this->rank(); 516 Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs)); 517 for (Index k = 0; k < rank; ++k) { 518 if (k != rank - 1) { 519 rhs.row(k).swap(rhs.row(rank - 1)); 520 } 521 rhs.middleRows(rank - 1, cols - rank + 1) 522 .applyHouseholderOnTheLeft( 523 matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k), 524 &temp(0)); 525 if (k != rank - 1) { 526 rhs.row(k).swap(rhs.row(rank - 1)); 527 } 528 } 529 } 530 531 #ifndef EIGEN_PARSED_BY_DOXYGEN 532 template <typename _MatrixType> 533 template <typename RhsType, typename DstType> 534 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl( 535 const RhsType& rhs, DstType& dst) const { 536 const Index rank = this->rank(); 537 if (rank == 0) { 538 dst.setZero(); 539 return; 540 } 541 542 // Compute c = Q^* * rhs 543 typename RhsType::PlainObject c(rhs); 544 c.applyOnTheLeft(matrixQ().setLength(rank).adjoint()); 545 546 // Solve T z = c(1:rank, :) 547 dst.topRows(rank) = matrixT() 548 .topLeftCorner(rank, rank) 549 .template triangularView<Upper>() 550 .solve(c.topRows(rank)); 551 552 const Index cols = this->cols(); 553 if (rank < cols) { 554 // Compute y = Z^* * [ z ] 555 // [ 0 ] 556 dst.bottomRows(cols - rank).setZero(); 557 applyZAdjointOnTheLeftInPlace(dst); 558 } 559 560 // Undo permutation to get x = P^{-1} * y. 561 dst = colsPermutation() * dst; 562 } 563 564 template<typename _MatrixType> 565 template<bool Conjugate, typename RhsType, typename DstType> 566 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const 567 { 568 const Index rank = this->rank(); 569 570 if (rank == 0) { 571 dst.setZero(); 572 return; 573 } 574 575 typename RhsType::PlainObject c(colsPermutation().transpose()*rhs); 576 577 if (rank < cols()) { 578 applyZOnTheLeftInPlace<!Conjugate>(c); 579 } 580 581 matrixT().topLeftCorner(rank, rank) 582 .template triangularView<Upper>() 583 .transpose().template conjugateIf<Conjugate>() 584 .solveInPlace(c.topRows(rank)); 585 586 dst.topRows(rank) = c.topRows(rank); 587 dst.bottomRows(rows()-rank).setZero(); 588 589 dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() ); 590 } 591 #endif 592 593 namespace internal { 594 595 template<typename MatrixType> 596 struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > > 597 : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject> 598 { 599 enum { Flags = 0 }; 600 }; 601 602 template<typename DstXprType, typename MatrixType> 603 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense> 604 { 605 typedef CompleteOrthogonalDecomposition<MatrixType> CodType; 606 typedef Inverse<CodType> SrcXprType; 607 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &) 608 { 609 typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType; 610 dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols())); 611 } 612 }; 613 614 } // end namespace internal 615 616 /** \returns the matrix Q as a sequence of householder transformations */ 617 template <typename MatrixType> 618 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType 619 CompleteOrthogonalDecomposition<MatrixType>::householderQ() const { 620 return m_cpqr.householderQ(); 621 } 622 623 /** \return the complete orthogonal decomposition of \c *this. 624 * 625 * \sa class CompleteOrthogonalDecomposition 626 */ 627 template <typename Derived> 628 const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject> 629 MatrixBase<Derived>::completeOrthogonalDecomposition() const { 630 return CompleteOrthogonalDecomposition<PlainObject>(eval()); 631 } 632 633 } // end namespace Eigen 634 635 #endif // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H