src/QR/FullPivHouseholderQR.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_FULLPIVOTINGHOUSEHOLDERQR_H
0012 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
0013
0014 namespace Eigen {
0015
0016 namespace internal {
0017
0018 template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
0019  : traits<_MatrixType>
0020 {
0021   typedef MatrixXpr XprKind;
0022   typedef SolverStorage StorageKind;
0023   typedef int StorageIndex;
0024   enum { Flags = 0 };
0025 };
0026
0027 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
0028
0029 template<typename MatrixType>
0030 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
0031 {
0032   typedef typename MatrixType::PlainObject ReturnType;
0033 };
0034
0035 } // end namespace internal
0036
0037 /** \ingroup QR_Module
0038   *
0039   * \class FullPivHouseholderQR
0040   *
0041   * \brief Householder rank-revealing QR decomposition of a matrix with full pivoting
0042   *
0043   * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
0044   *
0045   * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b P', \b Q and \b R
0046   * such that
0047   * \f[
0048   *  \mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R}
0049   * \f]
0050   * by using Householder transformations. Here, \b P and \b P' are permutation matrices, \b Q a unitary matrix
0051   * and \b R an upper triangular matrix.
0052   *
0053   * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal
0054   * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR.
0055   *
0056   * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
0057   *
0058   * \sa MatrixBase::fullPivHouseholderQr()
0059   */
0060 template<typename _MatrixType> class FullPivHouseholderQR
0061         : public SolverBase<FullPivHouseholderQR<_MatrixType> >
0062 {
0063   public:
0064
0065     typedef _MatrixType MatrixType;
0066     typedef SolverBase<FullPivHouseholderQR> Base;
0067     friend class SolverBase<FullPivHouseholderQR>;
0068
0069     EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivHouseholderQR)
0070     enum {
0071       MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
0072       MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
0073     };
0074     typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
0075     typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
0076     typedef Matrix<StorageIndex, 1,
0077                    EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
0078                    EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
0079     typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
0080     typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
0081     typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
0082     typedef typename MatrixType::PlainObject PlainObject;
0083
0084     /** \brief Default Constructor.
0085       *
0086       * The default constructor is useful in cases in which the user intends to
0087       * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
0088       */
0089     FullPivHouseholderQR()
0090       : m_qr(),
0091         m_hCoeffs(),
0092         m_rows_transpositions(),
0093         m_cols_transpositions(),
0094         m_cols_permutation(),
0095         m_temp(),
0096         m_isInitialized(false),
0097         m_usePrescribedThreshold(false) {}
0098
0099     /** \brief Default Constructor with memory preallocation
0100       *
0101       * Like the default constructor but with preallocation of the internal data
0102       * according to the specified problem \a size.
0103       * \sa FullPivHouseholderQR()
0104       */
0105     FullPivHouseholderQR(Index rows, Index cols)
0106       : m_qr(rows, cols),
0107         m_hCoeffs((std::min)(rows,cols)),
0108         m_rows_transpositions((std::min)(rows,cols)),
0109         m_cols_transpositions((std::min)(rows,cols)),
0110         m_cols_permutation(cols),
0111         m_temp(cols),
0112         m_isInitialized(false),
0113         m_usePrescribedThreshold(false) {}
0114
0115     /** \brief Constructs a QR factorization from a given matrix
0116       *
0117       * This constructor computes the QR factorization of the matrix \a matrix by calling
0118       * the method compute(). It is a short cut for:
0119       *
0120       * \code
0121       * FullPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
0122       * qr.compute(matrix);
0123       * \endcode
0124       *
0125       * \sa compute()
0126       */
0127     template<typename InputType>
0128     explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix)
0129       : m_qr(matrix.rows(), matrix.cols()),
0130         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
0131         m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
0132         m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
0133         m_cols_permutation(matrix.cols()),
0134         m_temp(matrix.cols()),
0135         m_isInitialized(false),
0136         m_usePrescribedThreshold(false)
0137     {
0138       compute(matrix.derived());
0139     }
0140
0141     /** \brief Constructs a QR factorization from a given matrix
0142       *
0143       * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
0144       *
0145       * \sa FullPivHouseholderQR(const EigenBase&)
0146       */
0147     template<typename InputType>
0148     explicit FullPivHouseholderQR(EigenBase<InputType>& matrix)
0149       : m_qr(matrix.derived()),
0150         m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
0151         m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
0152         m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
0153         m_cols_permutation(matrix.cols()),
0154         m_temp(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       * \c *this is the QR decomposition.
0164       *
0165       * \param b the right-hand-side of the equation to solve.
0166       *
0167       * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A,
0168       * and an arbitrary solution otherwise.
0169       *
0170       * \note_about_checking_solutions
0171       *
0172       * \note_about_arbitrary_choice_of_solution
0173       *
0174       * Example: \include FullPivHouseholderQR_solve.cpp
0175       * Output: \verbinclude FullPivHouseholderQR_solve.out
0176       */
0177     template<typename Rhs>
0178     inline const Solve<FullPivHouseholderQR, Rhs>
0179     solve(const MatrixBase<Rhs>& b) const;
0180     #endif
0181
0182     /** \returns Expression object representing the matrix Q
0183       */
0184     MatrixQReturnType matrixQ(void) const;
0185
0186     /** \returns a reference to the matrix where the Householder QR decomposition is stored
0187       */
0188     const MatrixType& matrixQR() const
0189     {
0190       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0191       return m_qr;
0192     }
0193
0194     template<typename InputType>
0195     FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
0196
0197     /** \returns a const reference to the column permutation matrix */
0198     const PermutationType& colsPermutation() const
0199     {
0200       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0201       return m_cols_permutation;
0202     }
0203
0204     /** \returns a const reference to the vector of indices representing the rows transpositions */
0205     const IntDiagSizeVectorType& rowsTranspositions() const
0206     {
0207       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0208       return m_rows_transpositions;
0209     }
0210
0211     /** \returns the absolute value of the determinant of the matrix of which
0212       * *this is the QR decomposition. It has only linear complexity
0213       * (that is, O(n) where n is the dimension of the square matrix)
0214       * as the QR decomposition has already been computed.
0215       *
0216       * \note This is only for square matrices.
0217       *
0218       * \warning a determinant can be very big or small, so for matrices
0219       * of large enough dimension, there is a risk of overflow/underflow.
0220       * One way to work around that is to use logAbsDeterminant() instead.
0221       *
0222       * \sa logAbsDeterminant(), MatrixBase::determinant()
0223       */
0224     typename MatrixType::RealScalar absDeterminant() const;
0225
0226     /** \returns the natural log of the absolute value of the determinant of the matrix of which
0227       * *this is the QR decomposition. It has only linear complexity
0228       * (that is, O(n) where n is the dimension of the square matrix)
0229       * as the QR decomposition has already been computed.
0230       *
0231       * \note This is only for square matrices.
0232       *
0233       * \note This method is useful to work around the risk of overflow/underflow that's inherent
0234       * to determinant computation.
0235       *
0236       * \sa absDeterminant(), MatrixBase::determinant()
0237       */
0238     typename MatrixType::RealScalar logAbsDeterminant() const;
0239
0240     /** \returns the rank of the matrix of which *this is the QR decomposition.
0241       *
0242       * \note This method has to determine which pivots should be considered nonzero.
0243       *       For that, it uses the threshold value that you can control by calling
0244       *       setThreshold(const RealScalar&).
0245       */
0246     inline Index rank() const
0247     {
0248       using std::abs;
0249       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0250       RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
0251       Index result = 0;
0252       for(Index i = 0; i < m_nonzero_pivots; ++i)
0253         result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
0254       return result;
0255     }
0256
0257     /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
0258       *
0259       * \note This method has to determine which pivots should be considered nonzero.
0260       *       For that, it uses the threshold value that you can control by calling
0261       *       setThreshold(const RealScalar&).
0262       */
0263     inline Index dimensionOfKernel() const
0264     {
0265       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0266       return cols() - rank();
0267     }
0268
0269     /** \returns true if the matrix of which *this is the QR decomposition represents an injective
0270       *          linear map, i.e. has trivial kernel; false otherwise.
0271       *
0272       * \note This method has to determine which pivots should be considered nonzero.
0273       *       For that, it uses the threshold value that you can control by calling
0274       *       setThreshold(const RealScalar&).
0275       */
0276     inline bool isInjective() const
0277     {
0278       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0279       return rank() == cols();
0280     }
0281
0282     /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
0283       *          linear map; false otherwise.
0284       *
0285       * \note This method has to determine which pivots should be considered nonzero.
0286       *       For that, it uses the threshold value that you can control by calling
0287       *       setThreshold(const RealScalar&).
0288       */
0289     inline bool isSurjective() const
0290     {
0291       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0292       return rank() == rows();
0293     }
0294
0295     /** \returns true if the matrix of which *this is the QR decomposition is invertible.
0296       *
0297       * \note This method has to determine which pivots should be considered nonzero.
0298       *       For that, it uses the threshold value that you can control by calling
0299       *       setThreshold(const RealScalar&).
0300       */
0301     inline bool isInvertible() const
0302     {
0303       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0304       return isInjective() && isSurjective();
0305     }
0306
0307     /** \returns the inverse of the matrix of which *this is the QR decomposition.
0308       *
0309       * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
0310       *       Use isInvertible() to first determine whether this matrix is invertible.
0311       */
0312     inline const Inverse<FullPivHouseholderQR> inverse() const
0313     {
0314       eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0315       return Inverse<FullPivHouseholderQR>(*this);
0316     }
0317
0318     inline Index rows() const { return m_qr.rows(); }
0319     inline Index cols() const { return m_qr.cols(); }
0320
0321     /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
0322       *
0323       * For advanced uses only.
0324       */
0325     const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
0326
0327     /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
0328       * who need to determine when pivots are to be considered nonzero. This is not used for the
0329       * QR decomposition itself.
0330       *
0331       * When it needs to get the threshold value, Eigen calls threshold(). By default, this
0332       * uses a formula to automatically determine a reasonable threshold.
0333       * Once you have called the present method setThreshold(const RealScalar&),
0334       * your value is used instead.
0335       *
0336       * \param threshold The new value to use as the threshold.
0337       *
0338       * A pivot will be considered nonzero if its absolute value is strictly greater than
0339       *  \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
0340       * where maxpivot is the biggest pivot.
0341       *
0342       * If you want to come back to the default behavior, call setThreshold(Default_t)
0343       */
0344     FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
0345     {
0346       m_usePrescribedThreshold = true;
0347       m_prescribedThreshold = threshold;
0348       return *this;
0349     }
0350
0351     /** Allows to come back to the default behavior, letting Eigen use its default formula for
0352       * determining the threshold.
0353       *
0354       * You should pass the special object Eigen::Default as parameter here.
0355       * \code qr.setThreshold(Eigen::Default); \endcode
0356       *
0357       * See the documentation of setThreshold(const RealScalar&).
0358       */
0359     FullPivHouseholderQR& setThreshold(Default_t)
0360     {
0361       m_usePrescribedThreshold = false;
0362       return *this;
0363     }
0364
0365     /** Returns the threshold that will be used by certain methods such as rank().
0366       *
0367       * See the documentation of setThreshold(const RealScalar&).
0368       */
0369     RealScalar threshold() const
0370     {
0371       eigen_assert(m_isInitialized || m_usePrescribedThreshold);
0372       return m_usePrescribedThreshold ? m_prescribedThreshold
0373       // this formula comes from experimenting (see "LU precision tuning" thread on the list)
0374       // and turns out to be identical to Higham's formula used already in LDLt.
0375                                       : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
0376     }
0377
0378     /** \returns the number of nonzero pivots in the QR decomposition.
0379       * Here nonzero is meant in the exact sense, not in a fuzzy sense.
0380       * So that notion isn't really intrinsically interesting, but it is
0381       * still useful when implementing algorithms.
0382       *
0383       * \sa rank()
0384       */
0385     inline Index nonzeroPivots() const
0386     {
0387       eigen_assert(m_isInitialized && "LU is not initialized.");
0388       return m_nonzero_pivots;
0389     }
0390
0391     /** \returns the absolute value of the biggest pivot, i.e. the biggest
0392       *          diagonal coefficient of U.
0393       */
0394     RealScalar maxPivot() const { return m_maxpivot; }
0395
0396     #ifndef EIGEN_PARSED_BY_DOXYGEN
0397     template<typename RhsType, typename DstType>
0398     void _solve_impl(const RhsType &rhs, DstType &dst) const;
0399
0400     template<bool Conjugate, typename RhsType, typename DstType>
0401     void _solve_impl_transposed(const RhsType &rhs, DstType &dst) const;
0402     #endif
0403
0404   protected:
0405
0406     static void check_template_parameters()
0407     {
0408       EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
0409     }
0410
0411     void computeInPlace();
0412
0413     MatrixType m_qr;
0414     HCoeffsType m_hCoeffs;
0415     IntDiagSizeVectorType m_rows_transpositions;
0416     IntDiagSizeVectorType m_cols_transpositions;
0417     PermutationType m_cols_permutation;
0418     RowVectorType m_temp;
0419     bool m_isInitialized, m_usePrescribedThreshold;
0420     RealScalar m_prescribedThreshold, m_maxpivot;
0421     Index m_nonzero_pivots;
0422     RealScalar m_precision;
0423     Index m_det_pq;
0424 };
0425
0426 template<typename MatrixType>
0427 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
0428 {
0429   using std::abs;
0430   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0431   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0432   return abs(m_qr.diagonal().prod());
0433 }
0434
0435 template<typename MatrixType>
0436 typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
0437 {
0438   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0439   eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
0440   return m_qr.diagonal().cwiseAbs().array().log().sum();
0441 }
0442
0443 /** Performs the QR factorization of the given matrix \a matrix. The result of
0444   * the factorization is stored into \c *this, and a reference to \c *this
0445   * is returned.
0446   *
0447   * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
0448   */
0449 template<typename MatrixType>
0450 template<typename InputType>
0451 FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
0452 {
0453   m_qr = matrix.derived();
0454   computeInPlace();
0455   return *this;
0456 }
0457
0458 template<typename MatrixType>
0459 void FullPivHouseholderQR<MatrixType>::computeInPlace()
0460 {
0461   check_template_parameters();
0462
0463   using std::abs;
0464   Index rows = m_qr.rows();
0465   Index cols = m_qr.cols();
0466   Index size = (std::min)(rows,cols);
0467
0468
0469   m_hCoeffs.resize(size);
0470
0471   m_temp.resize(cols);
0472
0473   m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
0474
0475   m_rows_transpositions.resize(size);
0476   m_cols_transpositions.resize(size);
0477   Index number_of_transpositions = 0;
0478
0479   RealScalar biggest(0);
0480
0481   m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
0482   m_maxpivot = RealScalar(0);
0483
0484   for (Index k = 0; k < size; ++k)
0485   {
0486     Index row_of_biggest_in_corner, col_of_biggest_in_corner;
0487     typedef internal::scalar_score_coeff_op<Scalar> Scoring;
0488     typedef typename Scoring::result_type Score;
0489
0490     Score score = m_qr.bottomRightCorner(rows-k, cols-k)
0491                       .unaryExpr(Scoring())
0492                       .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
0493     row_of_biggest_in_corner += k;
0494     col_of_biggest_in_corner += k;
0495     RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
0496     if(k==0) biggest = biggest_in_corner;
0497
0498     // if the corner is negligible, then we have less than full rank, and we can finish early
0499     if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
0500     {
0501       m_nonzero_pivots = k;
0502       for(Index i = k; i < size; i++)
0503       {
0504         m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
0505         m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
0506         m_hCoeffs.coeffRef(i) = Scalar(0);
0507       }
0508       break;
0509     }
0510
0511     m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner);
0512     m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner);
0513     if(k != row_of_biggest_in_corner) {
0514       m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
0515       ++number_of_transpositions;
0516     }
0517     if(k != col_of_biggest_in_corner) {
0518       m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
0519       ++number_of_transpositions;
0520     }
0521
0522     RealScalar beta;
0523     m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
0524     m_qr.coeffRef(k,k) = beta;
0525
0526     // remember the maximum absolute value of diagonal coefficients
0527     if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
0528
0529     m_qr.bottomRightCorner(rows-k, cols-k-1)
0530         .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
0531   }
0532
0533   m_cols_permutation.setIdentity(cols);
0534   for(Index k = 0; k < size; ++k)
0535     m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
0536
0537   m_det_pq = (number_of_transpositions%2) ? -1 : 1;
0538   m_isInitialized = true;
0539 }
0540
0541 #ifndef EIGEN_PARSED_BY_DOXYGEN
0542 template<typename _MatrixType>
0543 template<typename RhsType, typename DstType>
0544 void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
0545 {
0546   const Index l_rank = rank();
0547
0548   // FIXME introduce nonzeroPivots() and use it here. and more generally,
0549   // make the same improvements in this dec as in FullPivLU.
0550   if(l_rank==0)
0551   {
0552     dst.setZero();
0553     return;
0554   }
0555
0556   typename RhsType::PlainObject c(rhs);
0557
0558   Matrix<typename RhsType::Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
0559   for (Index k = 0; k < l_rank; ++k)
0560   {
0561     Index remainingSize = rows()-k;
0562     c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
0563     c.bottomRightCorner(remainingSize, rhs.cols())
0564       .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
0565                                m_hCoeffs.coeff(k), &temp.coeffRef(0));
0566   }
0567
0568   m_qr.topLeftCorner(l_rank, l_rank)
0569       .template triangularView<Upper>()
0570       .solveInPlace(c.topRows(l_rank));
0571
0572   for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
0573   for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
0574 }
0575
0576 template<typename _MatrixType>
0577 template<bool Conjugate, typename RhsType, typename DstType>
0578 void FullPivHouseholderQR<_MatrixType>::_solve_impl_transposed(const RhsType &rhs, DstType &dst) const
0579 {
0580   const Index l_rank = rank();
0581
0582   if(l_rank == 0)
0583   {
0584     dst.setZero();
0585     return;
0586   }
0587
0588   typename RhsType::PlainObject c(m_cols_permutation.transpose()*rhs);
0589
0590   m_qr.topLeftCorner(l_rank, l_rank)
0591          .template triangularView<Upper>()
0592          .transpose().template conjugateIf<Conjugate>()
0593          .solveInPlace(c.topRows(l_rank));
0594
0595   dst.topRows(l_rank) = c.topRows(l_rank);
0596   dst.bottomRows(rows()-l_rank).setZero();
0597
0598   Matrix<Scalar, 1, DstType::ColsAtCompileTime> temp(dst.cols());
0599   const Index size = (std::min)(rows(), cols());
0600   for (Index k = size-1; k >= 0; --k)
0601   {
0602     Index remainingSize = rows()-k;
0603
0604     dst.bottomRightCorner(remainingSize, dst.cols())
0605        .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1).template conjugateIf<!Conjugate>(),
0606                                   m_hCoeffs.template conjugateIf<Conjugate>().coeff(k), &temp.coeffRef(0));
0607
0608     dst.row(k).swap(dst.row(m_rows_transpositions.coeff(k)));
0609   }
0610 }
0611 #endif
0612
0613 namespace internal {
0614
0615 template<typename DstXprType, typename MatrixType>
0616 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
0617 {
0618   typedef FullPivHouseholderQR<MatrixType> QrType;
0619   typedef Inverse<QrType> SrcXprType;
0620   static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
0621   {
0622     dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
0623   }
0624 };
0625
0626 /** \ingroup QR_Module
0627   *
0628   * \brief Expression type for return value of FullPivHouseholderQR::matrixQ()
0629   *
0630   * \tparam MatrixType type of underlying dense matrix
0631   */
0632 template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
0633   : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
0634 {
0635 public:
0636   typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
0637   typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
0638   typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
0639                  MatrixType::MaxRowsAtCompileTime> WorkVectorType;
0640
0641   FullPivHouseholderQRMatrixQReturnType(const MatrixType&       qr,
0642                                         const HCoeffsType&      hCoeffs,
0643                                         const IntDiagSizeVectorType& rowsTranspositions)
0644     : m_qr(qr),
0645       m_hCoeffs(hCoeffs),
0646       m_rowsTranspositions(rowsTranspositions)
0647   {}
0648
0649   template <typename ResultType>
0650   void evalTo(ResultType& result) const
0651   {
0652     const Index rows = m_qr.rows();
0653     WorkVectorType workspace(rows);
0654     evalTo(result, workspace);
0655   }
0656
0657   template <typename ResultType>
0658   void evalTo(ResultType& result, WorkVectorType& workspace) const
0659   {
0660     using numext::conj;
0661     // compute the product H'_0 H'_1 ... H'_n-1,
0662     // where H_k is the k-th Householder transformation I - h_k v_k v_k'
0663     // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
0664     const Index rows = m_qr.rows();
0665     const Index cols = m_qr.cols();
0666     const Index size = (std::min)(rows, cols);
0667     workspace.resize(rows);
0668     result.setIdentity(rows, rows);
0669     for (Index k = size-1; k >= 0; k--)
0670     {
0671       result.block(k, k, rows-k, rows-k)
0672             .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
0673       result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
0674     }
0675   }
0676
0677   Index rows() const { return m_qr.rows(); }
0678   Index cols() const { return m_qr.rows(); }
0679
0680 protected:
0681   typename MatrixType::Nested m_qr;
0682   typename HCoeffsType::Nested m_hCoeffs;
0683   typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
0684 };
0685
0686 // template<typename MatrixType>
0687 // struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
0688 //  : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
0689 // {};
0690
0691 } // end namespace internal
0692
0693 template<typename MatrixType>
0694 inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
0695 {
0696   eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
0697   return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
0698 }
0699
0700 /** \return the full-pivoting Householder QR decomposition of \c *this.
0701   *
0702   * \sa class FullPivHouseholderQR
0703   */
0704 template<typename Derived>
0705 const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
0706 MatrixBase<Derived>::fullPivHouseholderQr() const
0707 {
0708   return FullPivHouseholderQR<PlainObject>(eval());
0709 }
0710
0711 } // end namespace Eigen
0712
0713 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H