src/QR/CompleteOrthogonalDecomposition.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2016 Rasmus Munk Larsen <rmlarsen@google.com>
0005 //
0006 // This Source Code Form is subject to the terms of the Mozilla
0007 // Public License v. 2.0. If a copy of the MPL was not distributed
0008 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0009
0010 #ifndef EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
0011 #define EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H
0012
0013 namespace Eigen {
0014
0015 namespace internal {
0016 template <typename _MatrixType>
0017 struct traits<CompleteOrthogonalDecomposition<_MatrixType> >
0018     : traits<_MatrixType> {
0019   typedef MatrixXpr XprKind;
0020   typedef SolverStorage StorageKind;
0021   typedef int StorageIndex;
0022   enum { Flags = 0 };
0023 };
0024
0025 }  // end namespace internal
0026
0027 /** \ingroup QR_Module
0028   *
0029   * \class CompleteOrthogonalDecomposition
0030   *
0031   * \brief Complete orthogonal decomposition (COD) of a matrix.
0032   *
0033   * \param MatrixType the type of the matrix of which we are computing the COD.
0034   *
0035   * This class performs a rank-revealing complete orthogonal decomposition of a
0036   * matrix  \b A into matrices \b P, \b Q, \b T, and \b Z such that
0037   * \f[
0038   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \,
0039   *                     \begin{bmatrix} \mathbf{T} &  \mathbf{0} \\
0040   *                                     \mathbf{0} & \mathbf{0} \end{bmatrix} \, \mathbf{Z}
0041   * \f]
0042   * by using Householder transformations. Here, \b P is a permutation matrix,
0043   * \b Q and \b Z are unitary matrices and \b T an upper triangular matrix of
0044   * size rank-by-rank. \b A may be rank deficient.
0045   *
0046   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0047   *
0048   * \sa MatrixBase::completeOrthogonalDecomposition()
0049   */
0050 template <typename _MatrixType> class CompleteOrthogonalDecomposition
0051           : public SolverBase<CompleteOrthogonalDecomposition<_MatrixType> >
0052 {
0053  public:
0054   typedef _MatrixType MatrixType;
0055   typedef SolverBase<CompleteOrthogonalDecomposition> Base;
0056
0057   template<typename Derived>
0058   friend struct internal::solve_assertion;
0059
0060   EIGEN_GENERIC_PUBLIC_INTERFACE(CompleteOrthogonalDecomposition)
0061   enum {
0062     MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0063     MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0064   };
0065   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
0066   typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>
0067       PermutationType;
0068   typedef typename internal::plain_row_type<MatrixType, Index>::type
0069       IntRowVectorType;
0070   typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
0071   typedef typename internal::plain_row_type<MatrixType, RealScalar>::type
0072       RealRowVectorType;
0073   typedef HouseholderSequence<
0074       MatrixType, typename internal::remove_all<
0075                       typename HCoeffsType::ConjugateReturnType>::type>
0076       HouseholderSequenceType;
0077   typedef typename MatrixType::PlainObject PlainObject;
0078
0079  private:
0080   typedef typename PermutationType::Index PermIndexType;
0081
0082  public:
0083   /**
0084    * \brief Default Constructor.
0085    *
0086    * The default constructor is useful in cases in which the user intends to
0087    * perform decompositions via
0088    * \c CompleteOrthogonalDecomposition::compute(const* MatrixType&).
0089    */
0090   CompleteOrthogonalDecomposition() : m_cpqr(), m_zCoeffs(), m_temp() {}
0091
0092   /** \brief Default Constructor with memory preallocation
0093    *
0094    * Like the default constructor but with preallocation of the internal data
0095    * according to the specified problem \a size.
0096    * \sa CompleteOrthogonalDecomposition()
0097    */
0098   CompleteOrthogonalDecomposition(Index rows, Index cols)
0099       : m_cpqr(rows, cols), m_zCoeffs((std::min)(rows, cols)), m_temp(cols) {}
0100
0101   /** \brief Constructs a complete orthogonal decomposition from a given
0102    * matrix.
0103    *
0104    * This constructor computes the complete orthogonal decomposition of the
0105    * matrix \a matrix by calling the method compute(). The default
0106    * threshold for rank determination will be used. It is a short cut for:
0107    *
0108    * \code
0109    * CompleteOrthogonalDecomposition<MatrixType> cod(matrix.rows(),
0110    *                                                 matrix.cols());
0111    * cod.setThreshold(Default);
0112    * cod.compute(matrix);
0113    * \endcode
0114    *
0115    * \sa compute()
0116    */
0117   template <typename InputType>
0118   explicit CompleteOrthogonalDecomposition(const EigenBase<InputType>& matrix)
0119       : m_cpqr(matrix.rows(), matrix.cols()),
0120         m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
0121         m_temp(matrix.cols())
0122   {
0123     compute(matrix.derived());
0124   }
0125
0126   /** \brief Constructs a complete orthogonal decomposition from a given matrix
0127     *
0128     * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
0129     *
0130     * \sa CompleteOrthogonalDecomposition(const EigenBase&)
0131     */
0132   template<typename InputType>
0133   explicit CompleteOrthogonalDecomposition(EigenBase<InputType>& matrix)
0134     : m_cpqr(matrix.derived()),
0135       m_zCoeffs((std::min)(matrix.rows(), matrix.cols())),
0136       m_temp(matrix.cols())
0137   {
0138     computeInPlace();
0139   }
0140
0141   #ifdef EIGEN_PARSED_BY_DOXYGEN
0142   /** This method computes the minimum-norm solution X to a least squares
0143    * problem \f[\mathrm{minimize} \|A X - B\|, \f] where \b A is the matrix of
0144    * which \c *this is the complete orthogonal decomposition.
0145    *
0146    * \param b the right-hand sides of the problem to solve.
0147    *
0148    * \returns a solution.
0149    *
0150    */
0151   template <typename Rhs>
0152   inline const Solve<CompleteOrthogonalDecomposition, Rhs> solve(
0153       const MatrixBase<Rhs>& b) const;
0154   #endif
0155
0156   HouseholderSequenceType householderQ(void) const;
0157   HouseholderSequenceType matrixQ(void) const { return m_cpqr.householderQ(); }
0158
0159   /** \returns the matrix \b Z.
0160    */
0161   MatrixType matrixZ() const {
0162     MatrixType Z = MatrixType::Identity(m_cpqr.cols(), m_cpqr.cols());
0163     applyZOnTheLeftInPlace<false>(Z);
0164     return Z;
0165   }
0166
0167   /** \returns a reference to the matrix where the complete orthogonal
0168    * decomposition is stored
0169    */
0170   const MatrixType& matrixQTZ() const { return m_cpqr.matrixQR(); }
0171
0172   /** \returns a reference to the matrix where the complete orthogonal
0173    * decomposition is stored.
0174    * \warning The strict lower part and \code cols() - rank() \endcode right
0175    * columns of this matrix contains internal values.
0176    * Only the upper triangular part should be referenced. To get it, use
0177    * \code matrixT().template triangularView<Upper>() \endcode
0178    * For rank-deficient matrices, use
0179    * \code
0180    * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
0181    * \endcode
0182    */
0183   const MatrixType& matrixT() const { return m_cpqr.matrixQR(); }
0184
0185   template <typename InputType>
0186   CompleteOrthogonalDecomposition& compute(const EigenBase<InputType>& matrix) {
0187     // Compute the column pivoted QR factorization A P = Q R.
0188     m_cpqr.compute(matrix);
0189     computeInPlace();
0190     return *this;
0191   }
0192
0193   /** \returns a const reference to the column permutation matrix */
0194   const PermutationType& colsPermutation() const {
0195     return m_cpqr.colsPermutation();
0196   }
0197
0198   /** \returns the absolute value of the determinant of the matrix of which
0199    * *this is the complete orthogonal decomposition. It has only linear
0200    * complexity (that is, O(n) where n is the dimension of the square matrix)
0201    * as the complete orthogonal decomposition has already been computed.
0202    *
0203    * \note This is only for square matrices.
0204    *
0205    * \warning a determinant can be very big or small, so for matrices
0206    * of large enough dimension, there is a risk of overflow/underflow.
0207    * One way to work around that is to use logAbsDeterminant() instead.
0208    *
0209    * \sa logAbsDeterminant(), MatrixBase::determinant()
0210    */
0211   typename MatrixType::RealScalar absDeterminant() const;
0212
0213   /** \returns the natural log of the absolute value of the determinant of the
0214    * matrix of which *this is the complete orthogonal decomposition. It has
0215    * only linear complexity (that is, O(n) where n is the dimension of the
0216    * square matrix) as the complete orthogonal decomposition has already been
0217    * computed.
0218    *
0219    * \note This is only for square matrices.
0220    *
0221    * \note This method is useful to work around the risk of overflow/underflow
0222    * that's inherent to determinant computation.
0223    *
0224    * \sa absDeterminant(), MatrixBase::determinant()
0225    */
0226   typename MatrixType::RealScalar logAbsDeterminant() const;
0227
0228   /** \returns the rank of the matrix of which *this is the complete orthogonal
0229    * decomposition.
0230    *
0231    * \note This method has to determine which pivots should be considered
0232    * nonzero. For that, it uses the threshold value that you can control by
0233    * calling setThreshold(const RealScalar&).
0234    */
0235   inline Index rank() const { return m_cpqr.rank(); }
0236
0237   /** \returns the dimension of the kernel of the matrix of which *this is the
0238    * complete orthogonal decomposition.
0239    *
0240    * \note This method has to determine which pivots should be considered
0241    * nonzero. For that, it uses the threshold value that you can control by
0242    * calling setThreshold(const RealScalar&).
0243    */
0244   inline Index dimensionOfKernel() const { return m_cpqr.dimensionOfKernel(); }
0245
0246   /** \returns true if the matrix of which *this is the decomposition represents
0247    * an injective linear map, i.e. has trivial kernel; false otherwise.
0248    *
0249    * \note This method has to determine which pivots should be considered
0250    * nonzero. For that, it uses the threshold value that you can control by
0251    * calling setThreshold(const RealScalar&).
0252    */
0253   inline bool isInjective() const { return m_cpqr.isInjective(); }
0254
0255   /** \returns true if the matrix of which *this is the decomposition represents
0256    * a surjective linear map; false otherwise.
0257    *
0258    * \note This method has to determine which pivots should be considered
0259    * nonzero. For that, it uses the threshold value that you can control by
0260    * calling setThreshold(const RealScalar&).
0261    */
0262   inline bool isSurjective() const { return m_cpqr.isSurjective(); }
0263
0264   /** \returns true if the matrix of which *this is the complete orthogonal
0265    * decomposition is invertible.
0266    *
0267    * \note This method has to determine which pivots should be considered
0268    * nonzero. For that, it uses the threshold value that you can control by
0269    * calling setThreshold(const RealScalar&).
0270    */
0271   inline bool isInvertible() const { return m_cpqr.isInvertible(); }
0272
0273   /** \returns the pseudo-inverse of the matrix of which *this is the complete
0274    * orthogonal decomposition.
0275    * \warning: Do not compute \c this->pseudoInverse()*rhs to solve a linear systems.
0276    * It is more efficient and numerically stable to call \c this->solve(rhs).
0277    */
0278   inline const Inverse<CompleteOrthogonalDecomposition> pseudoInverse() const
0279   {
0280     eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized.");
0281     return Inverse<CompleteOrthogonalDecomposition>(*this);
0282   }
0283
0284   inline Index rows() const { return m_cpqr.rows(); }
0285   inline Index cols() const { return m_cpqr.cols(); }
0286
0287   /** \returns a const reference to the vector of Householder coefficients used
0288    * to represent the factor \c Q.
0289    *
0290    * For advanced uses only.
0291    */
0292   inline const HCoeffsType& hCoeffs() const { return m_cpqr.hCoeffs(); }
0293
0294   /** \returns a const reference to the vector of Householder coefficients
0295    * used to represent the factor \c Z.
0296    *
0297    * For advanced uses only.
0298    */
0299   const HCoeffsType& zCoeffs() const { return m_zCoeffs; }
0300
0301   /** Allows to prescribe a threshold to be used by certain methods, such as
0302    * rank(), who need to determine when pivots are to be considered nonzero.
0303    * Most be called before calling compute().
0304    *
0305    * When it needs to get the threshold value, Eigen calls threshold(). By
0306    * default, this uses a formula to automatically determine a reasonable
0307    * threshold. Once you have called the present method
0308    * setThreshold(const RealScalar&), your value is used instead.
0309    *
0310    * \param threshold The new value to use as the threshold.
0311    *
0312    * A pivot will be considered nonzero if its absolute value is strictly
0313    * greater than
0314    *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
0315    * where maxpivot is the biggest pivot.
0316    *
0317    * If you want to come back to the default behavior, call
0318    * setThreshold(Default_t)
0319    */
0320   CompleteOrthogonalDecomposition& setThreshold(const RealScalar& threshold) {
0321     m_cpqr.setThreshold(threshold);
0322     return *this;
0323   }
0324
0325   /** Allows to come back to the default behavior, letting Eigen use its default
0326    * formula for determining the threshold.
0327    *
0328    * You should pass the special object Eigen::Default as parameter here.
0329    * \code qr.setThreshold(Eigen::Default); \endcode
0330    *
0331    * See the documentation of setThreshold(const RealScalar&).
0332    */
0333   CompleteOrthogonalDecomposition& setThreshold(Default_t) {
0334     m_cpqr.setThreshold(Default);
0335     return *this;
0336   }
0337
0338   /** Returns the threshold that will be used by certain methods such as rank().
0339    *
0340    * See the documentation of setThreshold(const RealScalar&).
0341    */
0342   RealScalar threshold() const { return m_cpqr.threshold(); }
0343
0344   /** \returns the number of nonzero pivots in the complete orthogonal
0345    * decomposition. Here nonzero is meant in the exact sense, not in a
0346    * fuzzy sense. So that notion isn't really intrinsically interesting,
0347    * but it is still useful when implementing algorithms.
0348    *
0349    * \sa rank()
0350    */
0351   inline Index nonzeroPivots() const { return m_cpqr.nonzeroPivots(); }
0352
0353   /** \returns the absolute value of the biggest pivot, i.e. the biggest
0354    *          diagonal coefficient of R.
0355    */
0356   inline RealScalar maxPivot() const { return m_cpqr.maxPivot(); }
0357
0358   /** \brief Reports whether the complete orthogonal decomposition was
0359    * successful.
0360    *
0361    * \note This function always returns \c Success. It is provided for
0362    * compatibility
0363    * with other factorization routines.
0364    * \returns \c Success
0365    */
0366   ComputationInfo info() const {
0367     eigen_assert(m_cpqr.m_isInitialized && "Decomposition is not initialized.");
0368     return Success;
0369   }
0370
0371 #ifndef EIGEN_PARSED_BY_DOXYGEN
0372   template <typename RhsType, typename DstType>
0373   void _solve_impl(const RhsType& rhs, DstType& dst) const;
0374
0375   template<bool Conjugate, typename RhsType, typename DstType>
0376   void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
0377 #endif
0378
0379  protected:
0380   static void check_template_parameters() {
0381     EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0382   }
0383
0384   template<bool Transpose_, typename Rhs>
0385   void _check_solve_assertion(const Rhs& b) const {
0386       EIGEN_ONLY_USED_FOR_DEBUG(b);
0387       eigen_assert(m_cpqr.m_isInitialized && "CompleteOrthogonalDecomposition is not initialized.");
0388       eigen_assert((Transpose_?derived().cols():derived().rows())==b.rows() && "CompleteOrthogonalDecomposition::solve(): invalid number of rows of the right hand side matrix b");
0389   }
0390
0391   void computeInPlace();
0392
0393   /** Overwrites \b rhs with \f$ \mathbf{Z} * \mathbf{rhs} \f$ or
0394    *  \f$ \mathbf{\overline Z} * \mathbf{rhs} \f$ if \c Conjugate
0395    *  is set to \c true.
0396    */
0397   template <bool Conjugate, typename Rhs>
0398   void applyZOnTheLeftInPlace(Rhs& rhs) const;
0399
0400   /** Overwrites \b rhs with \f$ \mathbf{Z}^* * \mathbf{rhs} \f$.
0401    */
0402   template <typename Rhs>
0403   void applyZAdjointOnTheLeftInPlace(Rhs& rhs) const;
0404
0405   ColPivHouseholderQR<MatrixType> m_cpqr;
0406   HCoeffsType m_zCoeffs;
0407   RowVectorType m_temp;
0408 };
0409
0410 template <typename MatrixType>
0411 typename MatrixType::RealScalar
0412 CompleteOrthogonalDecomposition<MatrixType>::absDeterminant() const {
0413   return m_cpqr.absDeterminant();
0414 }
0415
0416 template <typename MatrixType>
0417 typename MatrixType::RealScalar
0418 CompleteOrthogonalDecomposition<MatrixType>::logAbsDeterminant() const {
0419   return m_cpqr.logAbsDeterminant();
0420 }
0421
0422 /** Performs the complete orthogonal decomposition of the given matrix \a
0423  * matrix. The result of the factorization is stored into \c *this, and a
0424  * reference to \c *this is returned.
0425  *
0426  * \sa class CompleteOrthogonalDecomposition,
0427  * CompleteOrthogonalDecomposition(const MatrixType&)
0428  */
0429 template <typename MatrixType>
0430 void CompleteOrthogonalDecomposition<MatrixType>::computeInPlace()
0431 {
0432   check_template_parameters();
0433
0434   // the column permutation is stored as int indices, so just to be sure:
0435   eigen_assert(m_cpqr.cols() <= NumTraits<int>::highest());
0436
0437   const Index rank = m_cpqr.rank();
0438   const Index cols = m_cpqr.cols();
0439   const Index rows = m_cpqr.rows();
0440   m_zCoeffs.resize((std::min)(rows, cols));
0441   m_temp.resize(cols);
0442
0443   if (rank < cols) {
0444     // We have reduced the (permuted) matrix to the form
0445     //   [R11 R12]
0446     //   [ 0  R22]
0447     // where R11 is r-by-r (r = rank) upper triangular, R12 is
0448     // r-by-(n-r), and R22 is empty or the norm of R22 is negligible.
0449     // We now compute the complete orthogonal decomposition by applying
0450     // Householder transformations from the right to the upper trapezoidal
0451     // matrix X = [R11 R12] to zero out R12 and obtain the factorization
0452     // [R11 R12] = [T11 0] * Z, where T11 is r-by-r upper triangular and
0453     // Z = Z(0) * Z(1) ... Z(r-1) is an n-by-n orthogonal matrix.
0454     // We store the data representing Z in R12 and m_zCoeffs.
0455     for (Index k = rank - 1; k >= 0; --k) {
0456       if (k != rank - 1) {
0457         // Given the API for Householder reflectors, it is more convenient if
0458         // we swap the leading parts of columns k and r-1 (zero-based) to form
0459         // the matrix X_k = [X(0:k, k), X(0:k, r:n)]
0460         m_cpqr.m_qr.col(k).head(k + 1).swap(
0461             m_cpqr.m_qr.col(rank - 1).head(k + 1));
0462       }
0463       // Construct Householder reflector Z(k) to zero out the last row of X_k,
0464       // i.e. choose Z(k) such that
0465       // [X(k, k), X(k, r:n)] * Z(k) = [beta, 0, .., 0].
0466       RealScalar beta;
0467       m_cpqr.m_qr.row(k)
0468           .tail(cols - rank + 1)
0469           .makeHouseholderInPlace(m_zCoeffs(k), beta);
0470       m_cpqr.m_qr(k, rank - 1) = beta;
0471       if (k > 0) {
0472         // Apply Z(k) to the first k rows of X_k
0473         m_cpqr.m_qr.topRightCorner(k, cols - rank + 1)
0474             .applyHouseholderOnTheRight(
0475                 m_cpqr.m_qr.row(k).tail(cols - rank).adjoint(), m_zCoeffs(k),
0476                 &m_temp(0));
0477       }
0478       if (k != rank - 1) {
0479         // Swap X(0:k,k) back to its proper location.
0480         m_cpqr.m_qr.col(k).head(k + 1).swap(
0481             m_cpqr.m_qr.col(rank - 1).head(k + 1));
0482       }
0483     }
0484   }
0485 }
0486
0487 template <typename MatrixType>
0488 template <bool Conjugate, typename Rhs>
0489 void CompleteOrthogonalDecomposition<MatrixType>::applyZOnTheLeftInPlace(
0490     Rhs& rhs) const {
0491   const Index cols = this->cols();
0492   const Index nrhs = rhs.cols();
0493   const Index rank = this->rank();
0494   Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
0495   for (Index k = rank-1; k >= 0; --k) {
0496     if (k != rank - 1) {
0497       rhs.row(k).swap(rhs.row(rank - 1));
0498     }
0499     rhs.middleRows(rank - 1, cols - rank + 1)
0500         .applyHouseholderOnTheLeft(
0501             matrixQTZ().row(k).tail(cols - rank).transpose().template conjugateIf<!Conjugate>(), zCoeffs().template conjugateIf<Conjugate>()(k),
0502             &temp(0));
0503     if (k != rank - 1) {
0504       rhs.row(k).swap(rhs.row(rank - 1));
0505     }
0506   }
0507 }
0508
0509 template <typename MatrixType>
0510 template <typename Rhs>
0511 void CompleteOrthogonalDecomposition<MatrixType>::applyZAdjointOnTheLeftInPlace(
0512     Rhs& rhs) const {
0513   const Index cols = this->cols();
0514   const Index nrhs = rhs.cols();
0515   const Index rank = this->rank();
0516   Matrix<typename Rhs::Scalar, Dynamic, 1> temp((std::max)(cols, nrhs));
0517   for (Index k = 0; k < rank; ++k) {
0518     if (k != rank - 1) {
0519       rhs.row(k).swap(rhs.row(rank - 1));
0520     }
0521     rhs.middleRows(rank - 1, cols - rank + 1)
0522         .applyHouseholderOnTheLeft(
0523             matrixQTZ().row(k).tail(cols - rank).adjoint(), zCoeffs()(k),
0524             &temp(0));
0525     if (k != rank - 1) {
0526       rhs.row(k).swap(rhs.row(rank - 1));
0527     }
0528   }
0529 }
0530
0531 #ifndef EIGEN_PARSED_BY_DOXYGEN
0532 template <typename _MatrixType>
0533 template <typename RhsType, typename DstType>
0534 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl(
0535     const RhsType& rhs, DstType& dst) const {
0536   const Index rank = this->rank();
0537   if (rank == 0) {
0538     dst.setZero();
0539     return;
0540   }
0541
0542   // Compute c = Q^* * rhs
0543   typename RhsType::PlainObject c(rhs);
0544   c.applyOnTheLeft(matrixQ().setLength(rank).adjoint());
0545
0546   // Solve T z = c(1:rank, :)
0547   dst.topRows(rank) = matrixT()
0548                           .topLeftCorner(rank, rank)
0549                           .template triangularView<Upper>()
0550                           .solve(c.topRows(rank));
0551
0552   const Index cols = this->cols();
0553   if (rank < cols) {
0554     // Compute y = Z^* * [ z ]
0555     //                   [ 0 ]
0556     dst.bottomRows(cols - rank).setZero();
0557     applyZAdjointOnTheLeftInPlace(dst);
0558   }
0559
0560   // Undo permutation to get x = P^{-1} * y.
0561   dst = colsPermutation() * dst;
0562 }
0563
0564 template<typename _MatrixType>
0565 template<bool Conjugate, typename RhsType, typename DstType>
0566 void CompleteOrthogonalDecomposition<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
0567 {
0568   const Index rank = this->rank();
0569
0570   if (rank == 0) {
0571     dst.setZero();
0572     return;
0573   }
0574
0575   typename RhsType::PlainObject c(colsPermutation().transpose()*rhs);
0576
0577   if (rank < cols()) {
0578     applyZOnTheLeftInPlace<!Conjugate>(c);
0579   }
0580
0581   matrixT().topLeftCorner(rank, rank)
0582            .template triangularView<Upper>()
0583            .transpose().template conjugateIf<Conjugate>()
0584            .solveInPlace(c.topRows(rank));
0585
0586   dst.topRows(rank) = c.topRows(rank);
0587   dst.bottomRows(rows()-rank).setZero();
0588
0589   dst.applyOnTheLeft(householderQ().setLength(rank).template conjugateIf<!Conjugate>() );
0590 }
0591 #endif
0592
0593 namespace internal {
0594
0595 template<typename MatrixType>
0596 struct traits<Inverse<CompleteOrthogonalDecomposition<MatrixType> > >
0597   : traits<typename Transpose<typename MatrixType::PlainObject>::PlainObject>
0598 {
0599   enum { Flags = 0 };
0600 };
0601
0602 template<typename DstXprType, typename MatrixType>
0603 struct Assignment<DstXprType, Inverse<CompleteOrthogonalDecomposition<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename CompleteOrthogonalDecomposition<MatrixType>::Scalar>, Dense2Dense>
0604 {
0605   typedef CompleteOrthogonalDecomposition<MatrixType> CodType;
0606   typedef Inverse<CodType> SrcXprType;
0607   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename CodType::Scalar> &)
0608   {
0609     typedef Matrix<typename CodType::Scalar, CodType::RowsAtCompileTime, CodType::RowsAtCompileTime, 0, CodType::MaxRowsAtCompileTime, CodType::MaxRowsAtCompileTime> IdentityMatrixType;
0610     dst = src.nestedExpression().solve(IdentityMatrixType::Identity(src.cols(), src.cols()));
0611   }
0612 };
0613
0614 } // end namespace internal
0615
0616 /** \returns the matrix Q as a sequence of householder transformations */
0617 template <typename MatrixType>
0618 typename CompleteOrthogonalDecomposition<MatrixType>::HouseholderSequenceType
0619 CompleteOrthogonalDecomposition<MatrixType>::householderQ() const {
0620   return m_cpqr.householderQ();
0621 }
0622
0623 /** \return the complete orthogonal decomposition of \c *this.
0624   *
0625   * \sa class CompleteOrthogonalDecomposition
0626   */
0627 template <typename Derived>
0628 const CompleteOrthogonalDecomposition<typename MatrixBase<Derived>::PlainObject>
0629 MatrixBase<Derived>::completeOrthogonalDecomposition() const {
0630   return CompleteOrthogonalDecomposition<PlainObject>(eval());
0631 }
0632
0633 }  // end namespace Eigen
0634
0635 #endif  // EIGEN_COMPLETEORTHOGONALDECOMPOSITION_H