src/QR/ColPivHouseholderQR.h

0001 // This file is part of Eigen, a lightweight C++ template library
0002 // for linear algebra.
0003 //
0004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
0005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
0006 //
0007 // This Source Code Form is subject to the terms of the Mozilla
0008 // Public License v. 2.0. If a copy of the MPL was not distributed
0009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
0010
0011 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
0012 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017 template<typename _MatrixType> struct traits<ColPivHouseholderQR<_MatrixType> >
0018  : traits<_MatrixType>
0019 {
0020   typedef MatrixXpr XprKind;
0021   typedef SolverStorage StorageKind;
0022   typedef int StorageIndex;
0023   enum { Flags = 0 };
0024 };
0025
0026 } // end namespace internal
0027
0028 /** \ingroup QR_Module
0029   *
0030   * \class ColPivHouseholderQR
0031   *
0032   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
0033   *
0034   * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
0035   *
0036   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
0037   * such that
0038   * \f[
0039   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
0040   * \f]
0041   * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an
0042   * upper triangular matrix.
0043   *
0044   * This decomposition performs column pivoting in order to be rank-revealing and improve
0045   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
0046   *
0047   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0048   *
0049   * \sa MatrixBase::colPivHouseholderQr()
0050   */
0051 template<typename _MatrixType> class ColPivHouseholderQR
0052         : public SolverBase<ColPivHouseholderQR<_MatrixType> >
0053 {
0054   public:
0055
0056     typedef _MatrixType MatrixType;
0057     typedef SolverBase<ColPivHouseholderQR> Base;
0058     friend class SolverBase<ColPivHouseholderQR>;
0059
0060     EIGEN_GENERIC_PUBLIC_INTERFACE(ColPivHouseholderQR)
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> PermutationType;
0067     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
0068     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
0069     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
0070     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type> HouseholderSequenceType;
0071     typedef typename MatrixType::PlainObject PlainObject;
0072
0073   private:
0074
0075     typedef typename PermutationType::StorageIndex PermIndexType;
0076
0077   public:
0078
0079     /**
0080     * \brief Default Constructor.
0081     *
0082     * The default constructor is useful in cases in which the user intends to
0083     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
0084     */
0085     ColPivHouseholderQR()
0086       : m_qr(),
0087         m_hCoeffs(),
0088         m_colsPermutation(),
0089         m_colsTranspositions(),
0090         m_temp(),
0091         m_colNormsUpdated(),
0092         m_colNormsDirect(),
0093         m_isInitialized(false),
0094         m_usePrescribedThreshold(false) {}
0095
0096     /** \brief Default Constructor with memory preallocation
0097       *
0098       * Like the default constructor but with preallocation of the internal data
0099       * according to the specified problem \a size.
0100       * \sa ColPivHouseholderQR()
0101       */
0102     ColPivHouseholderQR(Index rows, Index cols)
0103       : m_qr(rows, cols),
0104         m_hCoeffs((std::min)(rows,cols)),
0105         m_colsPermutation(PermIndexType(cols)),
0106         m_colsTranspositions(cols),
0107         m_temp(cols),
0108         m_colNormsUpdated(cols),
0109         m_colNormsDirect(cols),
0110         m_isInitialized(false),
0111         m_usePrescribedThreshold(false) {}
0112
0113     /** \brief Constructs a QR factorization from a given matrix
0114       *
0115       * This constructor computes the QR factorization of the matrix \a matrix by calling
0116       * the method compute(). It is a short cut for:
0117       *
0118       * \code
0119       * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
0120       * qr.compute(matrix);
0121       * \endcode
0122       *
0123       * \sa compute()
0124       */
0125     template<typename InputType>
0126     explicit ColPivHouseholderQR(const EigenBase<InputType>& matrix)
0127       : m_qr(matrix.rows(), matrix.cols()),
0128         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
0129         m_colsPermutation(PermIndexType(matrix.cols())),
0130         m_colsTranspositions(matrix.cols()),
0131         m_temp(matrix.cols()),
0132         m_colNormsUpdated(matrix.cols()),
0133         m_colNormsDirect(matrix.cols()),
0134         m_isInitialized(false),
0135         m_usePrescribedThreshold(false)
0136     {
0137       compute(matrix.derived());
0138     }
0139
0140     /** \brief Constructs a QR factorization from a given matrix
0141       *
0142       * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
0143       *
0144       * \sa ColPivHouseholderQR(const EigenBase&)
0145       */
0146     template<typename InputType>
0147     explicit ColPivHouseholderQR(EigenBase<InputType>& matrix)
0148       : m_qr(matrix.derived()),
0149         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
0150         m_colsPermutation(PermIndexType(matrix.cols())),
0151         m_colsTranspositions(matrix.cols()),
0152         m_temp(matrix.cols()),
0153         m_colNormsUpdated(matrix.cols()),
0154         m_colNormsDirect(matrix.cols()),
0155         m_isInitialized(false),
0156         m_usePrescribedThreshold(false)
0157     {
0158       computeInPlace();
0159     }
0160
0161     #ifdef EIGEN_PARSED_BY_DOXYGEN
0162     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
0163       * *this is the QR decomposition, if any exists.
0164       *
0165       * \param b the right-hand-side of the equation to solve.
0166       *
0167       * \returns a solution.
0168       *
0169       * \note_about_checking_solutions
0170       *
0171       * \note_about_arbitrary_choice_of_solution
0172       *
0173       * Example: \include ColPivHouseholderQR_solve.cpp
0174       * Output: \verbinclude ColPivHouseholderQR_solve.out
0175       */
0176     template<typename Rhs>
0177     inline const Solve<ColPivHouseholderQR, Rhs>
0178     solve(const MatrixBase<Rhs>& b) const;
0179     #endif
0180
0181     HouseholderSequenceType householderQ() const;
0182     HouseholderSequenceType matrixQ() const
0183     {
0184       return householderQ();
0185     }
0186
0187     /** \returns a reference to the matrix where the Householder QR decomposition is stored
0188       */
0189     const MatrixType& matrixQR() const
0190     {
0191       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0192       return m_qr;
0193     }
0194
0195     /** \returns a reference to the matrix where the result Householder QR is stored
0196      * \warning The strict lower part of this matrix contains internal values.
0197      * Only the upper triangular part should be referenced. To get it, use
0198      * \code matrixR().template triangularView<Upper>() \endcode
0199      * For rank-deficient matrices, use
0200      * \code
0201      * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>()
0202      * \endcode
0203      */
0204     const MatrixType& matrixR() const
0205     {
0206       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0207       return m_qr;
0208     }
0209
0210     template<typename InputType>
0211     ColPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
0212
0213     /** \returns a const reference to the column permutation matrix */
0214     const PermutationType& colsPermutation() const
0215     {
0216       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0217       return m_colsPermutation;
0218     }
0219
0220     /** \returns the absolute value of the determinant of the matrix of which
0221       * *this is the QR decomposition. It has only linear complexity
0222       * (that is, O(n) where n is the dimension of the square matrix)
0223       * as the QR decomposition has already been computed.
0224       *
0225       * \note This is only for square matrices.
0226       *
0227       * \warning a determinant can be very big or small, so for matrices
0228       * of large enough dimension, there is a risk of overflow/underflow.
0229       * One way to work around that is to use logAbsDeterminant() instead.
0230       *
0231       * \sa logAbsDeterminant(), MatrixBase::determinant()
0232       */
0233     typename MatrixType::RealScalar absDeterminant() const;
0234
0235     /** \returns the natural log of the absolute value of the determinant of the matrix of which
0236       * *this is the QR decomposition. It has only linear complexity
0237       * (that is, O(n) where n is the dimension of the square matrix)
0238       * as the QR decomposition has already been computed.
0239       *
0240       * \note This is only for square matrices.
0241       *
0242       * \note This method is useful to work around the risk of overflow/underflow that's inherent
0243       * to determinant computation.
0244       *
0245       * \sa absDeterminant(), MatrixBase::determinant()
0246       */
0247     typename MatrixType::RealScalar logAbsDeterminant() const;
0248
0249     /** \returns the rank of the matrix of which *this is the QR decomposition.
0250       *
0251       * \note This method has to determine which pivots should be considered nonzero.
0252       *       For that, it uses the threshold value that you can control by calling
0253       *       setThreshold(const RealScalar&).
0254       */
0255     inline Index rank() const
0256     {
0257       using std::abs;
0258       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0259       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
0260       Index result = 0;
0261       for(Index i = 0; i < m_nonzero_pivots; ++i)
0262         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
0263       return result;
0264     }
0265
0266     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
0267       *
0268       * \note This method has to determine which pivots should be considered nonzero.
0269       *       For that, it uses the threshold value that you can control by calling
0270       *       setThreshold(const RealScalar&).
0271       */
0272     inline Index dimensionOfKernel() const
0273     {
0274       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0275       return cols() - rank();
0276     }
0277
0278     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
0279       *          linear map, i.e. has trivial kernel; false otherwise.
0280       *
0281       * \note This method has to determine which pivots should be considered nonzero.
0282       *       For that, it uses the threshold value that you can control by calling
0283       *       setThreshold(const RealScalar&).
0284       */
0285     inline bool isInjective() const
0286     {
0287       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0288       return rank() == cols();
0289     }
0290
0291     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
0292       *          linear map; false otherwise.
0293       *
0294       * \note This method has to determine which pivots should be considered nonzero.
0295       *       For that, it uses the threshold value that you can control by calling
0296       *       setThreshold(const RealScalar&).
0297       */
0298     inline bool isSurjective() const
0299     {
0300       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0301       return rank() == rows();
0302     }
0303
0304     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
0305       *
0306       * \note This method has to determine which pivots should be considered nonzero.
0307       *       For that, it uses the threshold value that you can control by calling
0308       *       setThreshold(const RealScalar&).
0309       */
0310     inline bool isInvertible() const
0311     {
0312       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0313       return isInjective() && isSurjective();
0314     }
0315
0316     /** \returns the inverse of the matrix of which *this is the QR decomposition.
0317       *
0318       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
0319       *       Use isInvertible() to first determine whether this matrix is invertible.
0320       */
0321     inline const Inverse<ColPivHouseholderQR> inverse() const
0322     {
0323       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0324       return Inverse<ColPivHouseholderQR>(*this);
0325     }
0326
0327     inline Index rows() const { return m_qr.rows(); }
0328     inline Index cols() const { return m_qr.cols(); }
0329
0330     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
0331       *
0332       * For advanced uses only.
0333       */
0334     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
0335
0336     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
0337       * who need to determine when pivots are to be considered nonzero. This is not used for the
0338       * QR decomposition itself.
0339       *
0340       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
0341       * uses a formula to automatically determine a reasonable threshold.
0342       * Once you have called the present method setThreshold(const RealScalar&),
0343       * your value is used instead.
0344       *
0345       * \param threshold The new value to use as the threshold.
0346       *
0347       * A pivot will be considered nonzero if its absolute value is strictly greater than
0348       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
0349       * where maxpivot is the biggest pivot.
0350       *
0351       * If you want to come back to the default behavior, call setThreshold(Default_t)
0352       */
0353     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
0354     {
0355       m_usePrescribedThreshold = true;
0356       m_prescribedThreshold = threshold;
0357       return *this;
0358     }
0359
0360     /** Allows to come back to the default behavior, letting Eigen use its default formula for
0361       * determining the threshold.
0362       *
0363       * You should pass the special object Eigen::Default as parameter here.
0364       * \code qr.setThreshold(Eigen::Default); \endcode
0365       *
0366       * See the documentation of setThreshold(const RealScalar&).
0367       */
0368     ColPivHouseholderQR& setThreshold(Default_t)
0369     {
0370       m_usePrescribedThreshold = false;
0371       return *this;
0372     }
0373
0374     /** Returns the threshold that will be used by certain methods such as rank().
0375       *
0376       * See the documentation of setThreshold(const RealScalar&).
0377       */
0378     RealScalar threshold() const
0379     {
0380       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
0381       return m_usePrescribedThreshold ? m_prescribedThreshold
0382       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
0383       // and turns out to be identical to Higham's formula used already in LDLt.
0384                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
0385     }
0386
0387     /** \returns the number of nonzero pivots in the QR decomposition.
0388       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
0389       * So that notion isn't really intrinsically interesting, but it is
0390       * still useful when implementing algorithms.
0391       *
0392       * \sa rank()
0393       */
0394     inline Index nonzeroPivots() const
0395     {
0396       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0397       return m_nonzero_pivots;
0398     }
0399
0400     /** \returns the absolute value of the biggest pivot, i.e. the biggest
0401       *          diagonal coefficient of R.
0402       */
0403     RealScalar maxPivot() const { return m_maxpivot; }
0404
0405     /** \brief Reports whether the QR factorization was successful.
0406       *
0407       * \note This function always returns \c Success. It is provided for compatibility
0408       * with other factorization routines.
0409       * \returns \c Success
0410       */
0411     ComputationInfo info() const
0412     {
0413       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
0414       return Success;
0415     }
0416
0417     #ifndef EIGEN_PARSED_BY_DOXYGEN
0418     template<typename RhsType, typename DstType>
0419     void _solve_impl(const RhsType &rhs, DstType &dst) const;
0420
0421     template<bool Conjugate, typename RhsType, typename DstType>
0422     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
0423     #endif
0424
0425   protected:
0426
0427     friend class CompleteOrthogonalDecomposition<MatrixType>;
0428
0429     static void check_template_parameters()
0430     {
0431       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0432     }
0433
0434     void computeInPlace();
0435
0436     MatrixType m_qr;
0437     HCoeffsType m_hCoeffs;
0438     PermutationType m_colsPermutation;
0439     IntRowVectorType m_colsTranspositions;
0440     RowVectorType m_temp;
0441     RealRowVectorType m_colNormsUpdated;
0442     RealRowVectorType m_colNormsDirect;
0443     bool m_isInitialized, m_usePrescribedThreshold;
0444     RealScalar m_prescribedThreshold, m_maxpivot;
0445     Index m_nonzero_pivots;
0446     Index m_det_pq;
0447 };
0448
0449 template<typename MatrixType>
0450 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant() const
0451 {
0452   using std::abs;
0453   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0454   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0455   return abs(m_qr.diagonal().prod());
0456 }
0457
0458 template<typename MatrixType>
0459 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
0460 {
0461   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0462   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0463   return m_qr.diagonal().cwiseAbs().array().log().sum();
0464 }
0465
0466 /** Performs the QR factorization of the given matrix \a matrix. The result of
0467   * the factorization is stored into \c *this, and a reference to \c *this
0468   * is returned.
0469   *
0470   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
0471   */
0472 template<typename MatrixType>
0473 template<typename InputType>
0474 ColPivHouseholderQR<MatrixType>& ColPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
0475 {
0476   m_qr = matrix.derived();
0477   computeInPlace();
0478   return *this;
0479 }
0480
0481 template<typename MatrixType>
0482 void ColPivHouseholderQR<MatrixType>::computeInPlace()
0483 {
0484   check_template_parameters();
0485
0486   // the column permutation is stored as int indices, so just to be sure:
0487   eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
0488
0489   using std::abs;
0490
0491   Index rows = m_qr.rows();
0492   Index cols = m_qr.cols();
0493   Index size = m_qr.diagonalSize();
0494
0495   m_hCoeffs.resize(size);
0496
0497   m_temp.resize(cols);
0498
0499   m_colsTranspositions.resize(m_qr.cols());
0500   Index number_of_transpositions = 0;
0501
0502   m_colNormsUpdated.resize(cols);
0503   m_colNormsDirect.resize(cols);
0504   for (Index k = 0; k < cols; ++k) {
0505     // colNormsDirect(k) caches the most recent directly computed norm of
0506     // column k.
0507     m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
0508     m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
0509   }
0510
0511   RealScalar threshold_helper =  numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
0512   RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
0513
0514   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
0515   m_maxpivot = RealScalar(0);
0516
0517   for(Index k = 0; k < size; ++k)
0518   {
0519     // first, we look up in our table m_colNormsUpdated which column has the biggest norm
0520     Index biggest_col_index;
0521     RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
0522     biggest_col_index += k;
0523
0524     // Track the number of meaningful pivots but do not stop the decomposition to make
0525     // sure that the initial matrix is properly reproduced. See bug 941.
0526     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
0527       m_nonzero_pivots = k;
0528
0529     // apply the transposition to the columns
0530     m_colsTranspositions.coeffRef(k) = biggest_col_index;
0531     if(k != biggest_col_index) {
0532       m_qr.col(k).swap(m_qr.col(biggest_col_index));
0533       std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
0534       std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
0535       ++number_of_transpositions;
0536     }
0537
0538     // generate the householder vector, store it below the diagonal
0539     RealScalar beta;
0540     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
0541
0542     // apply the householder transformation to the diagonal coefficient
0543     m_qr.coeffRef(k,k) = beta;
0544
0545     // remember the maximum absolute value of diagonal coefficients
0546     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
0547
0548     // apply the householder transformation
0549     m_qr.bottomRightCorner(rows-k, cols-k-1)
0550         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
0551
0552     // update our table of norms of the columns
0553     for (Index j = k + 1; j < cols; ++j) {
0554       // The following implements the stable norm downgrade step discussed in
0555       // http://www.netlib.org/lapack/lawnspdf/lawn176.pdf
0556       // and used in LAPACK routines xGEQPF and xGEQP3.
0557       // See lines 278-297 in http://www.netlib.org/lapack/explore-html/dc/df4/sgeqpf_8f_source.html
0558       if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
0559         RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
0560         temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
0561         temp = temp <  RealScalar(0) ? RealScalar(0) : temp;
0562         RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
0563                                                            m_colNormsDirect.coeffRef(j));
0564         if (temp2 <= norm_downdate_threshold) {
0565           // The updated norm has become too inaccurate so re-compute the column
0566           // norm directly.
0567           m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
0568           m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
0569         } else {
0570           m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
0571         }
0572       }
0573     }
0574   }
0575
0576   m_colsPermutation.setIdentity(PermIndexType(cols));
0577   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
0578     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
0579
0580   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
0581   m_isInitialized = true;
0582 }
0583
0584 #ifndef EIGEN_PARSED_BY_DOXYGEN
0585 template<typename _MatrixType>
0586 template<typename RhsType, typename DstType>
0587 void ColPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
0588 {
0589   const Index nonzero_pivots = nonzeroPivots();
0590
0591   if(nonzero_pivots == 0)
0592   {
0593     dst.setZero();
0594     return;
0595   }
0596
0597   typename RhsType::PlainObject c(rhs);
0598
0599   c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
0600
0601   m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
0602       .template triangularView<Upper>()
0603       .solveInPlace(c.topRows(nonzero_pivots));
0604
0605   for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
0606   for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
0607 }
0608
0609 template<typename _MatrixType>
0610 template<bool Conjugate, typename RhsType, typename DstType>
0611 void ColPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
0612 {
0613   const Index nonzero_pivots = nonzeroPivots();
0614
0615   if(nonzero_pivots == 0)
0616   {
0617     dst.setZero();
0618     return;
0619   }
0620
0621   typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs);
0622
0623   m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
0624         .template triangularView<Upper>()
0625         .transpose().template conjugateIf<Conjugate>()
0626         .solveInPlace(c.topRows(nonzero_pivots));
0627
0628   dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
0629   dst.bottomRows(rows()-nonzero_pivots).setZero();
0630
0631   dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).template conjugateIf<!Conjugate>() );
0632 }
0633 #endif
0634
0635 namespace internal {
0636
0637 template<typename DstXprType, typename MatrixType>
0638 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
0639 {
0640   typedef ColPivHouseholderQR<MatrixType> QrType;
0641   typedef Inverse<QrType> SrcXprType;
0642   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
0643   {
0644     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
0645   }
0646 };
0647
0648 } // end namespace internal
0649
0650 /** \returns the matrix Q as a sequence of householder transformations.
0651   * You can extract the meaningful part only by using:
0652   * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
0653 template<typename MatrixType>
0654 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
0655   ::householderQ() const
0656 {
0657   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
0658   return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
0659 }
0660
0661 /** \return the column-pivoting Householder QR decomposition of \c *this.
0662   *
0663   * \sa class ColPivHouseholderQR
0664   */
0665 template<typename Derived>
0666 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
0667 MatrixBase<Derived>::colPivHouseholderQr() const
0668 {
0669   return ColPivHouseholderQR<PlainObject>(eval());
0670 }
0671
0672 } // end namespace Eigen
0673
0674 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H