ColPivHouseholderQR.h

00001 // This file is part of Eigen, a lightweight C++ template library
00002 // for linear algebra.
00003 //
00004 // Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
00005 // Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
00006 //
00007 // This Source Code Form is subject to the terms of the Mozilla
00008 // Public License v. 2.0. If a copy of the MPL was not distributed
00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
00010 
00011 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00012 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
00013 
00014 namespace Eigen { 
00015 
00016 /** \ingroup QR_Module
00017   *
00018   * \class ColPivHouseholderQR
00019   *
00020   * \brief Householder rank-revealing QR decomposition of a matrix with column-pivoting
00021   *
00022   * \param MatrixType the type of the matrix of which we are computing the QR decomposition
00023   *
00024   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b Q and \b R
00025   * such that 
00026   * \f[
00027   *  \mathbf{A} \, \mathbf{P} = \mathbf{Q} \, \mathbf{R}
00028   * \f]
00029   * by using Householder transformations. Here, \b P is a permutation matrix, \b Q a unitary matrix and \b R an 
00030   * upper triangular matrix.
00031   *
00032   * This decomposition performs column pivoting in order to be rank-revealing and improve
00033   * numerical stability. It is slower than HouseholderQR, and faster than FullPivHouseholderQR.
00034   *
00035   * \sa MatrixBase::colPivHouseholderQr()
00036   */
00037 template<typename _MatrixType> class ColPivHouseholderQR
00038 {
00039   public:
00040 
00041     typedef _MatrixType MatrixType;
00042     enum {
00043       RowsAtCompileTime = MatrixType::RowsAtCompileTime,
00044       ColsAtCompileTime = MatrixType::ColsAtCompileTime,
00045       Options = MatrixType::Options,
00046       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
00047       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
00048     };
00049     typedef typename MatrixType::Scalar Scalar;
00050     typedef typename MatrixType::RealScalar RealScalar;
00051     typedef typename MatrixType::Index Index;
00052     typedef Matrix<Scalar, RowsAtCompileTime, RowsAtCompileTime, Options, MaxRowsAtCompileTime, MaxRowsAtCompileTime> MatrixQType;
00053     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
00054     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime>  PermutationType ;
00055     typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
00056     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
00057     typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
00058     typedef HouseholderSequence<MatrixType,typename internal::remove_all<typename HCoeffsType::ConjugateReturnType>::type > HouseholderSequenceType ;
00059     
00060   private:
00061     
00062     typedef typename PermutationType::Index PermIndexType;
00063     
00064   public:
00065 
00066     /**
00067     * \brief Default Constructor.
00068     *
00069     * The default constructor is useful in cases in which the user intends to
00070     * perform decompositions via ColPivHouseholderQR::compute(const MatrixType&).
00071     */
00072     ColPivHouseholderQR()
00073       : m_qr(),
00074         m_hCoeffs(),
00075         m_colsPermutation(),
00076         m_colsTranspositions(),
00077         m_temp(),
00078         m_colSqNorms(),
00079         m_isInitialized(false),
00080         m_usePrescribedThreshold(false) {}
00081 
00082     /** \brief Default Constructor with memory preallocation
00083       *
00084       * Like the default constructor but with preallocation of the internal data
00085       * according to the specified problem \a size.
00086       * \sa ColPivHouseholderQR()
00087       */
00088     ColPivHouseholderQR(Index rows, Index cols)
00089       : m_qr(rows, cols),
00090         m_hCoeffs((std::min)(rows,cols)),
00091         m_colsPermutation(PermIndexType(cols)),
00092         m_colsTranspositions(cols),
00093         m_temp(cols),
00094         m_colSqNorms(cols),
00095         m_isInitialized(false),
00096         m_usePrescribedThreshold(false) {}
00097 
00098     /** \brief Constructs a QR factorization from a given matrix
00099       *
00100       * This constructor computes the QR factorization of the matrix \a matrix by calling
00101       * the method compute(). It is a short cut for:
00102       * 
00103       * \code
00104       * ColPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
00105       * qr.compute(matrix);
00106       * \endcode
00107       * 
00108       * \sa compute()
00109       */
00110     ColPivHouseholderQR(const MatrixType& matrix)
00111       : m_qr(matrix.rows(), matrix.cols()),
00112         m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
00113         m_colsPermutation(PermIndexType(matrix.cols())),
00114         m_colsTranspositions(matrix.cols()),
00115         m_temp(matrix.cols()),
00116         m_colSqNorms(matrix.cols()),
00117         m_isInitialized(false),
00118         m_usePrescribedThreshold(false)
00119     {
00120       compute(matrix);
00121     }
00122 
00123     /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
00124       * *this is the QR decomposition, if any exists.
00125       *
00126       * \param b the right-hand-side of the equation to solve.
00127       *
00128       * \returns a solution.
00129       *
00130       * \note The case where b is a matrix is not yet implemented. Also, this
00131       *       code is space inefficient.
00132       *
00133       * \note_about_checking_solutions
00134       *
00135       * \note_about_arbitrary_choice_of_solution
00136       *
00137       * Example: \include ColPivHouseholderQR_solve.cpp
00138       * Output: \verbinclude ColPivHouseholderQR_solve.out
00139       */
00140     template<typename Rhs>
00141     inline const internal::solve_retval<ColPivHouseholderQR, Rhs>
00142     solve(const MatrixBase<Rhs>& b) const
00143     {
00144       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00145       return internal::solve_retval<ColPivHouseholderQR, Rhs>(*this, b.derived());
00146     }
00147 
00148     HouseholderSequenceType householderQ (void) const;
00149     HouseholderSequenceType matrixQ(void) const
00150     {
00151       return householderQ (); 
00152     }
00153 
00154     /** \returns a reference to the matrix where the Householder QR decomposition is stored
00155       */
00156     const MatrixType& matrixQR () const
00157     {
00158       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00159       return m_qr;
00160     }
00161     
00162     /** \returns a reference to the matrix where the result Householder QR is stored 
00163      * \warning The strict lower part of this matrix contains internal values. 
00164      * Only the upper triangular part should be referenced. To get it, use
00165      * \code matrixR().template triangularView<Upper>() \endcode
00166      * For rank-deficient matrices, use 
00167      * \code 
00168      * matrixR().topLeftCorner(rank(), rank()).template triangularView<Upper>() 
00169      * \endcode
00170      */
00171     const MatrixType& matrixR () const
00172     {
00173       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00174       return m_qr;
00175     }
00176     
00177     ColPivHouseholderQR& compute(const MatrixType& matrix);
00178 
00179     /** \returns a const reference to the column permutation matrix */
00180     const PermutationType & colsPermutation () const
00181     {
00182       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00183       return m_colsPermutation;
00184     }
00185 
00186     /** \returns the absolute value of the determinant of the matrix of which
00187       * *this is the QR decomposition. It has only linear complexity
00188       * (that is, O(n) where n is the dimension of the square matrix)
00189       * as the QR decomposition has already been computed.
00190       *
00191       * \note This is only for square matrices.
00192       *
00193       * \warning a determinant can be very big or small, so for matrices
00194       * of large enough dimension, there is a risk of overflow/underflow.
00195       * One way to work around that is to use logAbsDeterminant() instead.
00196       *
00197       * \sa logAbsDeterminant(), MatrixBase::determinant()
00198       */
00199     typename MatrixType::RealScalar absDeterminant () const;
00200 
00201     /** \returns the natural log of the absolute value of the determinant of the matrix of which
00202       * *this is the QR decomposition. It has only linear complexity
00203       * (that is, O(n) where n is the dimension of the square matrix)
00204       * as the QR decomposition has already been computed.
00205       *
00206       * \note This is only for square matrices.
00207       *
00208       * \note This method is useful to work around the risk of overflow/underflow that's inherent
00209       * to determinant computation.
00210       *
00211       * \sa absDeterminant(), MatrixBase::determinant()
00212       */
00213     typename MatrixType::RealScalar logAbsDeterminant () const;
00214 
00215     /** \returns the rank of the matrix of which *this is the QR decomposition.
00216       *
00217       * \note This method has to determine which pivots should be considered nonzero.
00218       *       For that, it uses the threshold value that you can control by calling
00219       *       setThreshold(const RealScalar&).
00220       */
00221     inline Index rank () const
00222     {
00223       using std::abs;
00224       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00225       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
00226       Index result = 0;
00227       for(Index i = 0; i < m_nonzero_pivots; ++i)
00228         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
00229       return result;
00230     }
00231 
00232     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
00233       *
00234       * \note This method has to determine which pivots should be considered nonzero.
00235       *       For that, it uses the threshold value that you can control by calling
00236       *       setThreshold(const RealScalar&).
00237       */
00238     inline Index dimensionOfKernel () const
00239     {
00240       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00241       return cols() - rank ();
00242     }
00243 
00244     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
00245       *          linear map, i.e. has trivial kernel; false otherwise.
00246       *
00247       * \note This method has to determine which pivots should be considered nonzero.
00248       *       For that, it uses the threshold value that you can control by calling
00249       *       setThreshold(const RealScalar&).
00250       */
00251     inline bool isInjective () const
00252     {
00253       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00254       return rank () == cols();
00255     }
00256 
00257     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
00258       *          linear map; false otherwise.
00259       *
00260       * \note This method has to determine which pivots should be considered nonzero.
00261       *       For that, it uses the threshold value that you can control by calling
00262       *       setThreshold(const RealScalar&).
00263       */
00264     inline bool isSurjective () const
00265     {
00266       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00267       return rank () == rows();
00268     }
00269 
00270     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
00271       *
00272       * \note This method has to determine which pivots should be considered nonzero.
00273       *       For that, it uses the threshold value that you can control by calling
00274       *       setThreshold(const RealScalar&).
00275       */
00276     inline bool isInvertible () const
00277     {
00278       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00279       return isInjective () && isSurjective ();
00280     }
00281 
00282     /** \returns the inverse of the matrix of which *this is the QR decomposition.
00283       *
00284       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
00285       *       Use isInvertible() to first determine whether this matrix is invertible.
00286       */
00287     inline const
00288     internal::solve_retval<ColPivHouseholderQR, typename MatrixType::IdentityReturnType>
00289     inverse () const
00290     {
00291       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00292       return internal::solve_retval<ColPivHouseholderQR,typename MatrixType::IdentityReturnType>
00293                (*this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
00294     }
00295 
00296     inline Index rows() const { return m_qr.rows(); }
00297     inline Index cols() const { return m_qr.cols(); }
00298     
00299     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
00300       * 
00301       * For advanced uses only.
00302       */
00303     const HCoeffsType& hCoeffs () const { return m_hCoeffs; }
00304 
00305     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
00306       * who need to determine when pivots are to be considered nonzero. This is not used for the
00307       * QR decomposition itself.
00308       *
00309       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
00310       * uses a formula to automatically determine a reasonable threshold.
00311       * Once you have called the present method setThreshold(const RealScalar&),
00312       * your value is used instead.
00313       *
00314       * \param threshold The new value to use as the threshold.
00315       *
00316       * A pivot will be considered nonzero if its absolute value is strictly greater than
00317       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
00318       * where maxpivot is the biggest pivot.
00319       *
00320       * If you want to come back to the default behavior, call setThreshold(Default_t)
00321       */
00322     ColPivHouseholderQR& setThreshold(const RealScalar& threshold)
00323     {
00324       m_usePrescribedThreshold = true;
00325       m_prescribedThreshold = threshold;
00326       return *this;
00327     }
00328 
00329     /** Allows to come back to the default behavior, letting Eigen use its default formula for
00330       * determining the threshold.
00331       *
00332       * You should pass the special object Eigen::Default as parameter here.
00333       * \code qr.setThreshold(Eigen::Default); \endcode
00334       *
00335       * See the documentation of setThreshold(const RealScalar&).
00336       */
00337     ColPivHouseholderQR& setThreshold(Default_t)
00338     {
00339       m_usePrescribedThreshold = false;
00340       return *this;
00341     }
00342 
00343     /** Returns the threshold that will be used by certain methods such as rank().
00344       *
00345       * See the documentation of setThreshold(const RealScalar&).
00346       */
00347     RealScalar threshold() const
00348     {
00349       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
00350       return m_usePrescribedThreshold ? m_prescribedThreshold
00351       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
00352       // and turns out to be identical to Higham's formula used already in LDLt.
00353                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
00354     }
00355 
00356     /** \returns the number of nonzero pivots in the QR decomposition.
00357       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
00358       * So that notion isn't really intrinsically interesting, but it is
00359       * still useful when implementing algorithms.
00360       *
00361       * \sa rank()
00362       */
00363     inline Index nonzeroPivots () const
00364     {
00365       eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00366       return m_nonzero_pivots;
00367     }
00368 
00369     /** \returns the absolute value of the biggest pivot, i.e. the biggest
00370       *          diagonal coefficient of R.
00371       */
00372     RealScalar maxPivot () const { return m_maxpivot; }
00373     
00374     /** \brief Reports whether the QR factorization was succesful.
00375       *
00376       * \note This function always returns \c Success. It is provided for compatibility 
00377       * with other factorization routines.
00378       * \returns \c Success 
00379       */
00380     ComputationInfo info() const
00381     {
00382       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00383       return Success;
00384     }
00385 
00386   protected:
00387     
00388     static void check_template_parameters()
00389     {
00390       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
00391     }
00392     
00393     MatrixType m_qr;
00394     HCoeffsType m_hCoeffs;
00395     PermutationType m_colsPermutation;
00396     IntRowVectorType m_colsTranspositions;
00397     RowVectorType m_temp;
00398     RealRowVectorType m_colSqNorms;
00399     bool m_isInitialized, m_usePrescribedThreshold;
00400     RealScalar m_prescribedThreshold, m_maxpivot;
00401     Index m_nonzero_pivots;
00402     Index m_det_pq;
00403 };
00404 
00405 template<typename MatrixType>
00406 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::absDeterminant () const
00407 {
00408   using std::abs;
00409   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00410   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00411   return abs(m_qr.diagonal().prod());
00412 }
00413 
00414 template<typename MatrixType>
00415 typename MatrixType::RealScalar ColPivHouseholderQR<MatrixType>::logAbsDeterminant() const
00416 {
00417   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00418   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
00419   return m_qr.diagonal().cwiseAbs().array().log().sum();
00420 }
00421 
00422 /** Performs the QR factorization of the given matrix \a matrix. The result of
00423   * the factorization is stored into \c *this, and a reference to \c *this
00424   * is returned.
00425   *
00426   * \sa class ColPivHouseholderQR, ColPivHouseholderQR(const MatrixType&)
00427   */
00428 template<typename MatrixType>
00429 ColPivHouseholderQR<MatrixType> & ColPivHouseholderQR<MatrixType>::compute(const MatrixType& matrix)
00430 {
00431   check_template_parameters();
00432   
00433   using std::abs;
00434   Index rows = matrix.rows();
00435   Index cols = matrix.cols();
00436   Index size = matrix.diagonalSize();
00437   
00438   // the column permutation is stored as int indices, so just to be sure:
00439   eigen_assert(cols<=NumTraits<int>::highest());
00440 
00441   m_qr = matrix;
00442   m_hCoeffs.resize(size);
00443 
00444   m_temp.resize(cols);
00445 
00446   m_colsTranspositions.resize(matrix.cols());
00447   Index number_of_transpositions = 0;
00448 
00449   m_colSqNorms.resize(cols);
00450   for(Index k = 0; k < cols; ++k)
00451     m_colSqNorms.coeffRef(k) = m_qr.col(k).squaredNorm();
00452 
00453   RealScalar threshold_helper = m_colSqNorms.maxCoeff() * numext::abs2(NumTraits<Scalar>::epsilon()) / RealScalar(rows);
00454 
00455   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
00456   m_maxpivot = RealScalar(0);
00457 
00458   for(Index k = 0; k < size; ++k)
00459   {
00460     // first, we look up in our table m_colSqNorms which column has the biggest squared norm
00461     Index biggest_col_index;
00462     RealScalar biggest_col_sq_norm = m_colSqNorms.tail(cols-k).maxCoeff(&biggest_col_index);
00463     biggest_col_index += k;
00464 
00465     // since our table m_colSqNorms accumulates imprecision at every step, we must now recompute
00466     // the actual squared norm of the selected column.
00467     // Note that not doing so does result in solve() sometimes returning inf/nan values
00468     // when running the unit test with 1000 repetitions.
00469     biggest_col_sq_norm = m_qr.col(biggest_col_index).tail(rows-k).squaredNorm();
00470 
00471     // we store that back into our table: it can't hurt to correct our table.
00472     m_colSqNorms.coeffRef(biggest_col_index) = biggest_col_sq_norm;
00473 
00474     // Track the number of meaningful pivots but do not stop the decomposition to make
00475     // sure that the initial matrix is properly reproduced. See bug 941.
00476     if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
00477       m_nonzero_pivots = k;
00478 
00479     // apply the transposition to the columns
00480     m_colsTranspositions.coeffRef(k) = biggest_col_index;
00481     if(k != biggest_col_index) {
00482       m_qr.col(k).swap(m_qr.col(biggest_col_index));
00483       std::swap(m_colSqNorms.coeffRef(k), m_colSqNorms.coeffRef(biggest_col_index));
00484       ++number_of_transpositions;
00485     }
00486 
00487     // generate the householder vector, store it below the diagonal
00488     RealScalar beta;
00489     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
00490 
00491     // apply the householder transformation to the diagonal coefficient
00492     m_qr.coeffRef(k,k) = beta;
00493 
00494     // remember the maximum absolute value of diagonal coefficients
00495     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
00496 
00497     // apply the householder transformation
00498     m_qr.bottomRightCorner(rows-k, cols-k-1)
00499         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
00500 
00501     // update our table of squared norms of the columns
00502     m_colSqNorms.tail(cols-k-1) -= m_qr.row(k).tail(cols-k-1).cwiseAbs2();
00503   }
00504 
00505   m_colsPermutation.setIdentity(PermIndexType(cols));
00506   for(PermIndexType k = 0; k < size/*m_nonzero_pivots*/; ++k)
00507     m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
00508 
00509   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
00510   m_isInitialized = true;
00511 
00512   return *this;
00513 }
00514 
00515 namespace internal {
00516 
00517 template<typename _MatrixType, typename Rhs>
00518 struct solve_retval<ColPivHouseholderQR<_MatrixType>, Rhs>
00519   : solve_retval_base<ColPivHouseholderQR<_MatrixType>, Rhs>
00520 {
00521   EIGEN_MAKE_SOLVE_HELPERS(ColPivHouseholderQR<_MatrixType>,Rhs)
00522 
00523   template<typename Dest> void evalTo(Dest& dst) const
00524   {
00525     eigen_assert(rhs().rows() == dec().rows());
00526 
00527     const Index cols = dec().cols(),
00528                 nonzero_pivots = dec().nonzeroPivots();
00529 
00530     if(nonzero_pivots == 0)
00531     {
00532       dst.setZero();
00533       return;
00534     }
00535 
00536     typename Rhs::PlainObject c(rhs());
00537 
00538     // Note that the matrix Q = H_0^* H_1^*... so its inverse is Q^* = (H_0 H_1 ...)^T
00539     c.applyOnTheLeft(householderSequence(dec().matrixQR(), dec().hCoeffs())
00540                      .setLength(dec().nonzeroPivots())
00541              .transpose()
00542       );
00543 
00544     dec().matrixR()
00545        .topLeftCorner(nonzero_pivots, nonzero_pivots)
00546        .template triangularView<Upper>()
00547        .solveInPlace(c.topRows(nonzero_pivots));
00548 
00549     for(Index i = 0; i < nonzero_pivots; ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
00550     for(Index i = nonzero_pivots; i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
00551   }
00552 };
00553 
00554 } // end namespace internal
00555 
00556 /** \returns the matrix Q as a sequence of householder transformations.
00557   * You can extract the meaningful part only by using:
00558   * \code qr.householderQ().setLength(qr.nonzeroPivots()) \endcode*/
00559 template<typename MatrixType>
00560 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
00561   ::householderQ() const
00562 {
00563   eigen_assert(m_isInitialized && "ColPivHouseholderQR is not initialized.");
00564   return HouseholderSequenceType (m_qr, m_hCoeffs.conjugate());
00565 }
00566 
00567 /** \return the column-pivoting Householder QR decomposition of \c *this.
00568   *
00569   * \sa class ColPivHouseholderQR
00570   */
00571 template<typename Derived>
00572 const ColPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
00573 MatrixBase<Derived>::colPivHouseholderQr() const
00574 {
00575   return ColPivHouseholderQR<PlainObject>(eval());
00576 }
00577 
00578 } // end namespace Eigen
00579 
00580 #endif // EIGEN_COLPIVOTINGHOUSEHOLDERQR_H