1// This file is part of Eigen, a lightweight C++ template library
2// for linear algebra.
3//
4// Copyright (C) 2008-2009 Gael Guennebaud <gael.guennebaud@inria.fr>
5// Copyright (C) 2009 Benoit Jacob <jacob.benoit.1@gmail.com>
6//
7// This Source Code Form is subject to the terms of the Mozilla
8// Public License v. 2.0. If a copy of the MPL was not distributed
9// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11#ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
12#define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
13
14namespace Eigen {
15
16namespace internal {
17
18template<typename _MatrixType> struct traits<FullPivHouseholderQR<_MatrixType> >
19 : traits<_MatrixType>
20{
21 enum { Flags = 0 };
22};
23
24template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType;
25
26template<typename MatrixType>
27struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
28{
29 typedef typename MatrixType::PlainObject ReturnType;
30};
31
32} // end namespace internal
33
34/** \ingroup QR_Module
35 *
36 * \class FullPivHouseholderQR
37 *
38 * \brief Householder rank-revealing QR decomposition of a matrix with full pivoting
39 *
40 * \tparam _MatrixType the type of the matrix of which we are computing the QR decomposition
41 *
42 * This class performs a rank-revealing QR decomposition of a matrix \b A into matrices \b P, \b P', \b Q and \b R
43 * such that
44 * \f[
45 * \mathbf{P} \, \mathbf{A} \, \mathbf{P}' = \mathbf{Q} \, \mathbf{R}
46 * \f]
47 * by using Householder transformations. Here, \b P and \b P' are permutation matrices, \b Q a unitary matrix
48 * and \b R an upper triangular matrix.
49 *
50 * This decomposition performs a very prudent full pivoting in order to be rank-revealing and achieve optimal
51 * numerical stability. The trade-off is that it is slower than HouseholderQR and ColPivHouseholderQR.
52 *
53 * This class supports the \link InplaceDecomposition inplace decomposition \endlink mechanism.
54 *
55 * \sa MatrixBase::fullPivHouseholderQr()
56 */
57template<typename _MatrixType> class FullPivHouseholderQR
58{
59 public:
60
61 typedef _MatrixType MatrixType;
62 enum {
63 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
64 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
65 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
66 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67 };
68 typedef typename MatrixType::Scalar Scalar;
69 typedef typename MatrixType::RealScalar RealScalar;
70 // FIXME should be int
71 typedef typename MatrixType::StorageIndex StorageIndex;
72 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
73 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74 typedef Matrix<StorageIndex, 1,
75 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime), RowMajor, 1,
76 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
77 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
78 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
79 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
80 typedef typename MatrixType::PlainObject PlainObject;
81
82 /** \brief Default Constructor.
83 *
84 * The default constructor is useful in cases in which the user intends to
85 * perform decompositions via FullPivHouseholderQR::compute(const MatrixType&).
86 */
87 FullPivHouseholderQR()
88 : m_qr(),
89 m_hCoeffs(),
90 m_rows_transpositions(),
91 m_cols_transpositions(),
92 m_cols_permutation(),
93 m_temp(),
94 m_isInitialized(false),
95 m_usePrescribedThreshold(false) {}
96
97 /** \brief Default Constructor with memory preallocation
98 *
99 * Like the default constructor but with preallocation of the internal data
100 * according to the specified problem \a size.
101 * \sa FullPivHouseholderQR()
102 */
103 FullPivHouseholderQR(Index rows, Index cols)
104 : m_qr(rows, cols),
105 m_hCoeffs((std::min)(rows,cols)),
106 m_rows_transpositions((std::min)(rows,cols)),
107 m_cols_transpositions((std::min)(rows,cols)),
108 m_cols_permutation(cols),
109 m_temp(cols),
110 m_isInitialized(false),
111 m_usePrescribedThreshold(false) {}
112
113 /** \brief Constructs a QR factorization from a given matrix
114 *
115 * This constructor computes the QR factorization of the matrix \a matrix by calling
116 * the method compute(). It is a short cut for:
117 *
118 * \code
119 * FullPivHouseholderQR<MatrixType> qr(matrix.rows(), matrix.cols());
120 * qr.compute(matrix);
121 * \endcode
122 *
123 * \sa compute()
124 */
125 template<typename InputType>
126 explicit FullPivHouseholderQR(const EigenBase<InputType>& matrix)
127 : m_qr(matrix.rows(), matrix.cols()),
128 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
129 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
130 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
131 m_cols_permutation(matrix.cols()),
132 m_temp(matrix.cols()),
133 m_isInitialized(false),
134 m_usePrescribedThreshold(false)
135 {
136 compute(matrix.derived());
137 }
138
139 /** \brief Constructs a QR factorization from a given matrix
140 *
141 * This overloaded constructor is provided for \link InplaceDecomposition inplace decomposition \endlink when \c MatrixType is a Eigen::Ref.
142 *
143 * \sa FullPivHouseholderQR(const EigenBase&)
144 */
145 template<typename InputType>
146 explicit FullPivHouseholderQR(EigenBase<InputType>& matrix)
147 : m_qr(matrix.derived()),
148 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
149 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
150 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
151 m_cols_permutation(matrix.cols()),
152 m_temp(matrix.cols()),
153 m_isInitialized(false),
154 m_usePrescribedThreshold(false)
155 {
156 computeInPlace();
157 }
158
159 /** This method finds a solution x to the equation Ax=b, where A is the matrix of which
160 * \c *this is the QR decomposition.
161 *
162 * \param b the right-hand-side of the equation to solve.
163 *
164 * \returns the exact or least-square solution if the rank is greater or equal to the number of columns of A,
165 * and an arbitrary solution otherwise.
166 *
167 * \note_about_checking_solutions
168 *
169 * \note_about_arbitrary_choice_of_solution
170 *
171 * Example: \include FullPivHouseholderQR_solve.cpp
172 * Output: \verbinclude FullPivHouseholderQR_solve.out
173 */
174 template<typename Rhs>
175 inline const Solve<FullPivHouseholderQR, Rhs>
176 solve(const MatrixBase<Rhs>& b) const
177 {
178 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
179 return Solve<FullPivHouseholderQR, Rhs>(*this, b.derived());
180 }
181
182 /** \returns Expression object representing the matrix Q
183 */
184 MatrixQReturnType matrixQ(void) const;
185
186 /** \returns a reference to the matrix where the Householder QR decomposition is stored
187 */
188 const MatrixType& matrixQR() const
189 {
190 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
191 return m_qr;
192 }
193
194 template<typename InputType>
195 FullPivHouseholderQR& compute(const EigenBase<InputType>& matrix);
196
197 /** \returns a const reference to the column permutation matrix */
198 const PermutationType& colsPermutation() const
199 {
200 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
201 return m_cols_permutation;
202 }
203
204 /** \returns a const reference to the vector of indices representing the rows transpositions */
205 const IntDiagSizeVectorType& rowsTranspositions() const
206 {
207 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
208 return m_rows_transpositions;
209 }
210
211 /** \returns the absolute value of the determinant of the matrix of which
212 * *this is the QR decomposition. It has only linear complexity
213 * (that is, O(n) where n is the dimension of the square matrix)
214 * as the QR decomposition has already been computed.
215 *
216 * \note This is only for square matrices.
217 *
218 * \warning a determinant can be very big or small, so for matrices
219 * of large enough dimension, there is a risk of overflow/underflow.
220 * One way to work around that is to use logAbsDeterminant() instead.
221 *
222 * \sa logAbsDeterminant(), MatrixBase::determinant()
223 */
224 typename MatrixType::RealScalar absDeterminant() const;
225
226 /** \returns the natural log of the absolute value of the determinant of the matrix of which
227 * *this is the QR decomposition. It has only linear complexity
228 * (that is, O(n) where n is the dimension of the square matrix)
229 * as the QR decomposition has already been computed.
230 *
231 * \note This is only for square matrices.
232 *
233 * \note This method is useful to work around the risk of overflow/underflow that's inherent
234 * to determinant computation.
235 *
236 * \sa absDeterminant(), MatrixBase::determinant()
237 */
238 typename MatrixType::RealScalar logAbsDeterminant() const;
239
240 /** \returns the rank of the matrix of which *this is the QR decomposition.
241 *
242 * \note This method has to determine which pivots should be considered nonzero.
243 * For that, it uses the threshold value that you can control by calling
244 * setThreshold(const RealScalar&).
245 */
246 inline Index rank() const
247 {
248 using std::abs;
249 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
250 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
251 Index result = 0;
252 for(Index i = 0; i < m_nonzero_pivots; ++i)
253 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
254 return result;
255 }
256
257 /** \returns the dimension of the kernel of the matrix of which *this is the QR decomposition.
258 *
259 * \note This method has to determine which pivots should be considered nonzero.
260 * For that, it uses the threshold value that you can control by calling
261 * setThreshold(const RealScalar&).
262 */
263 inline Index dimensionOfKernel() const
264 {
265 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
266 return cols() - rank();
267 }
268
269 /** \returns true if the matrix of which *this is the QR decomposition represents an injective
270 * linear map, i.e. has trivial kernel; false otherwise.
271 *
272 * \note This method has to determine which pivots should be considered nonzero.
273 * For that, it uses the threshold value that you can control by calling
274 * setThreshold(const RealScalar&).
275 */
276 inline bool isInjective() const
277 {
278 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
279 return rank() == cols();
280 }
281
282 /** \returns true if the matrix of which *this is the QR decomposition represents a surjective
283 * linear map; false otherwise.
284 *
285 * \note This method has to determine which pivots should be considered nonzero.
286 * For that, it uses the threshold value that you can control by calling
287 * setThreshold(const RealScalar&).
288 */
289 inline bool isSurjective() const
290 {
291 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
292 return rank() == rows();
293 }
294
295 /** \returns true if the matrix of which *this is the QR decomposition is invertible.
296 *
297 * \note This method has to determine which pivots should be considered nonzero.
298 * For that, it uses the threshold value that you can control by calling
299 * setThreshold(const RealScalar&).
300 */
301 inline bool isInvertible() const
302 {
303 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
304 return isInjective() && isSurjective();
305 }
306
307 /** \returns the inverse of the matrix of which *this is the QR decomposition.
308 *
309 * \note If this matrix is not invertible, the returned matrix has undefined coefficients.
310 * Use isInvertible() to first determine whether this matrix is invertible.
311 */
312 inline const Inverse<FullPivHouseholderQR> inverse() const
313 {
314 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
315 return Inverse<FullPivHouseholderQR>(*this);
316 }
317
318 inline Index rows() const { return m_qr.rows(); }
319 inline Index cols() const { return m_qr.cols(); }
320
321 /** \returns a const reference to the vector of Householder coefficients used to represent the factor \c Q.
322 *
323 * For advanced uses only.
324 */
325 const HCoeffsType& hCoeffs() const { return m_hCoeffs; }
326
327 /** Allows to prescribe a threshold to be used by certain methods, such as rank(),
328 * who need to determine when pivots are to be considered nonzero. This is not used for the
329 * QR decomposition itself.
330 *
331 * When it needs to get the threshold value, Eigen calls threshold(). By default, this
332 * uses a formula to automatically determine a reasonable threshold.
333 * Once you have called the present method setThreshold(const RealScalar&),
334 * your value is used instead.
335 *
336 * \param threshold The new value to use as the threshold.
337 *
338 * A pivot will be considered nonzero if its absolute value is strictly greater than
339 * \f$ \vert pivot \vert \leqslant threshold \times \vert maxpivot \vert \f$
340 * where maxpivot is the biggest pivot.
341 *
342 * If you want to come back to the default behavior, call setThreshold(Default_t)
343 */
344 FullPivHouseholderQR& setThreshold(const RealScalar& threshold)
345 {
346 m_usePrescribedThreshold = true;
347 m_prescribedThreshold = threshold;
348 return *this;
349 }
350
351 /** Allows to come back to the default behavior, letting Eigen use its default formula for
352 * determining the threshold.
353 *
354 * You should pass the special object Eigen::Default as parameter here.
355 * \code qr.setThreshold(Eigen::Default); \endcode
356 *
357 * See the documentation of setThreshold(const RealScalar&).
358 */
359 FullPivHouseholderQR& setThreshold(Default_t)
360 {
361 m_usePrescribedThreshold = false;
362 return *this;
363 }
364
365 /** Returns the threshold that will be used by certain methods such as rank().
366 *
367 * See the documentation of setThreshold(const RealScalar&).
368 */
369 RealScalar threshold() const
370 {
371 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
372 return m_usePrescribedThreshold ? m_prescribedThreshold
373 // this formula comes from experimenting (see "LU precision tuning" thread on the list)
374 // and turns out to be identical to Higham's formula used already in LDLt.
375 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
376 }
377
378 /** \returns the number of nonzero pivots in the QR decomposition.
379 * Here nonzero is meant in the exact sense, not in a fuzzy sense.
380 * So that notion isn't really intrinsically interesting, but it is
381 * still useful when implementing algorithms.
382 *
383 * \sa rank()
384 */
385 inline Index nonzeroPivots() const
386 {
387 eigen_assert(m_isInitialized && "LU is not initialized.");
388 return m_nonzero_pivots;
389 }
390
391 /** \returns the absolute value of the biggest pivot, i.e. the biggest
392 * diagonal coefficient of U.
393 */
394 RealScalar maxPivot() const { return m_maxpivot; }
395
396 #ifndef EIGEN_PARSED_BY_DOXYGEN
397 template<typename RhsType, typename DstType>
398 EIGEN_DEVICE_FUNC
399 void _solve_impl(const RhsType &rhs, DstType &dst) const;
400 #endif
401
402 protected:
403
404 static void check_template_parameters()
405 {
406 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
407 }
408
409 void computeInPlace();
410
411 MatrixType m_qr;
412 HCoeffsType m_hCoeffs;
413 IntDiagSizeVectorType m_rows_transpositions;
414 IntDiagSizeVectorType m_cols_transpositions;
415 PermutationType m_cols_permutation;
416 RowVectorType m_temp;
417 bool m_isInitialized, m_usePrescribedThreshold;
418 RealScalar m_prescribedThreshold, m_maxpivot;
419 Index m_nonzero_pivots;
420 RealScalar m_precision;
421 Index m_det_pq;
422};
423
424template<typename MatrixType>
425typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::absDeterminant() const
426{
427 using std::abs;
428 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
429 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
430 return abs(m_qr.diagonal().prod());
431}
432
433template<typename MatrixType>
434typename MatrixType::RealScalar FullPivHouseholderQR<MatrixType>::logAbsDeterminant() const
435{
436 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
437 eigen_assert(m_qr.rows() == m_qr.cols() && "You can't take the determinant of a non-square matrix!");
438 return m_qr.diagonal().cwiseAbs().array().log().sum();
439}
440
441/** Performs the QR factorization of the given matrix \a matrix. The result of
442 * the factorization is stored into \c *this, and a reference to \c *this
443 * is returned.
444 *
445 * \sa class FullPivHouseholderQR, FullPivHouseholderQR(const MatrixType&)
446 */
447template<typename MatrixType>
448template<typename InputType>
449FullPivHouseholderQR<MatrixType>& FullPivHouseholderQR<MatrixType>::compute(const EigenBase<InputType>& matrix)
450{
451 m_qr = matrix.derived();
452 computeInPlace();
453 return *this;
454}
455
456template<typename MatrixType>
457void FullPivHouseholderQR<MatrixType>::computeInPlace()
458{
459 check_template_parameters();
460
461 using std::abs;
462 Index rows = m_qr.rows();
463 Index cols = m_qr.cols();
464 Index size = (std::min)(rows,cols);
465
466
467 m_hCoeffs.resize(size);
468
469 m_temp.resize(cols);
470
471 m_precision = NumTraits<Scalar>::epsilon() * RealScalar(size);
472
473 m_rows_transpositions.resize(size);
474 m_cols_transpositions.resize(size);
475 Index number_of_transpositions = 0;
476
477 RealScalar biggest(0);
478
479 m_nonzero_pivots = size; // the generic case is that in which all pivots are nonzero (invertible case)
480 m_maxpivot = RealScalar(0);
481
482 for (Index k = 0; k < size; ++k)
483 {
484 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
485 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
486 typedef typename Scoring::result_type Score;
487
488 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
489 .unaryExpr(Scoring())
490 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
491 row_of_biggest_in_corner += k;
492 col_of_biggest_in_corner += k;
493 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
494 if(k==0) biggest = biggest_in_corner;
495
496 // if the corner is negligible, then we have less than full rank, and we can finish early
497 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
498 {
499 m_nonzero_pivots = k;
500 for(Index i = k; i < size; i++)
501 {
502 m_rows_transpositions.coeffRef(i) = i;
503 m_cols_transpositions.coeffRef(i) = i;
504 m_hCoeffs.coeffRef(i) = Scalar(0);
505 }
506 break;
507 }
508
509 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
510 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
511 if(k != row_of_biggest_in_corner) {
512 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
513 ++number_of_transpositions;
514 }
515 if(k != col_of_biggest_in_corner) {
516 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
517 ++number_of_transpositions;
518 }
519
520 RealScalar beta;
521 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
522 m_qr.coeffRef(k,k) = beta;
523
524 // remember the maximum absolute value of diagonal coefficients
525 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
526
527 m_qr.bottomRightCorner(rows-k, cols-k-1)
528 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
529 }
530
531 m_cols_permutation.setIdentity(cols);
532 for(Index k = 0; k < size; ++k)
533 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
534
535 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
536 m_isInitialized = true;
537}
538
539#ifndef EIGEN_PARSED_BY_DOXYGEN
540template<typename _MatrixType>
541template<typename RhsType, typename DstType>
542void FullPivHouseholderQR<_MatrixType>::_solve_impl(const RhsType &rhs, DstType &dst) const
543{
544 eigen_assert(rhs.rows() == rows());
545 const Index l_rank = rank();
546
547 // FIXME introduce nonzeroPivots() and use it here. and more generally,
548 // make the same improvements in this dec as in FullPivLU.
549 if(l_rank==0)
550 {
551 dst.setZero();
552 return;
553 }
554
555 typename RhsType::PlainObject c(rhs);
556
557 Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
558 for (Index k = 0; k < l_rank; ++k)
559 {
560 Index remainingSize = rows()-k;
561 c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
562 c.bottomRightCorner(remainingSize, rhs.cols())
563 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
564 m_hCoeffs.coeff(k), &temp.coeffRef(0));
565 }
566
567 m_qr.topLeftCorner(l_rank, l_rank)
568 .template triangularView<Upper>()
569 .solveInPlace(c.topRows(l_rank));
570
571 for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
572 for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
573}
574#endif
575
576namespace internal {
577
578template<typename DstXprType, typename MatrixType>
579struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
580{
581 typedef FullPivHouseholderQR<MatrixType> QrType;
582 typedef Inverse<QrType> SrcXprType;
583 static void run(DstXprType &dst, const SrcXprType &src, const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
584 {
585 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
586 }
587};
588
589/** \ingroup QR_Module
590 *
591 * \brief Expression type for return value of FullPivHouseholderQR::matrixQ()
592 *
593 * \tparam MatrixType type of underlying dense matrix
594 */
595template<typename MatrixType> struct FullPivHouseholderQRMatrixQReturnType
596 : public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
597{
598public:
599 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
600 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
601 typedef Matrix<typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1,
602 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
603
604 FullPivHouseholderQRMatrixQReturnType(const MatrixType& qr,
605 const HCoeffsType& hCoeffs,
606 const IntDiagSizeVectorType& rowsTranspositions)
607 : m_qr(qr),
608 m_hCoeffs(hCoeffs),
609 m_rowsTranspositions(rowsTranspositions)
610 {}
611
612 template <typename ResultType>
613 void evalTo(ResultType& result) const
614 {
615 const Index rows = m_qr.rows();
616 WorkVectorType workspace(rows);
617 evalTo(result, workspace);
618 }
619
620 template <typename ResultType>
621 void evalTo(ResultType& result, WorkVectorType& workspace) const
622 {
623 using numext::conj;
624 // compute the product H'_0 H'_1 ... H'_n-1,
625 // where H_k is the k-th Householder transformation I - h_k v_k v_k'
626 // and v_k is the k-th Householder vector [1,m_qr(k+1,k), m_qr(k+2,k), ...]
627 const Index rows = m_qr.rows();
628 const Index cols = m_qr.cols();
629 const Index size = (std::min)(rows, cols);
630 workspace.resize(rows);
631 result.setIdentity(rows, rows);
632 for (Index k = size-1; k >= 0; k--)
633 {
634 result.block(k, k, rows-k, rows-k)
635 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
636 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
637 }
638 }
639
640 Index rows() const { return m_qr.rows(); }
641 Index cols() const { return m_qr.rows(); }
642
643protected:
644 typename MatrixType::Nested m_qr;
645 typename HCoeffsType::Nested m_hCoeffs;
646 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
647};
648
649// template<typename MatrixType>
650// struct evaluator<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
651// : public evaluator<ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> > >
652// {};
653
654} // end namespace internal
655
656template<typename MatrixType>
657inline typename FullPivHouseholderQR<MatrixType>::MatrixQReturnType FullPivHouseholderQR<MatrixType>::matrixQ() const
658{
659 eigen_assert(m_isInitialized && "FullPivHouseholderQR is not initialized.");
660 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
661}
662
663/** \return the full-pivoting Householder QR decomposition of \c *this.
664 *
665 * \sa class FullPivHouseholderQR
666 */
667template<typename Derived>
668const FullPivHouseholderQR<typename MatrixBase<Derived>::PlainObject>
669MatrixBase<Derived>::fullPivHouseholderQr() const
670{
671 return FullPivHouseholderQR<PlainObject>(eval());
672}
673
674} // end namespace Eigen
675
676#endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
677