11 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 18 template<
typename _MatrixType>
struct traits<FullPivHouseholderQR<_MatrixType> >
24 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType;
26 template<
typename MatrixType>
27 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
29 typedef typename MatrixType::PlainObject ReturnType;
55 template<
typename _MatrixType>
class FullPivHouseholderQR
59 typedef _MatrixType MatrixType;
61 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
62 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
63 Options = MatrixType::Options,
64 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67 typedef typename MatrixType::Scalar Scalar;
68 typedef typename MatrixType::RealScalar RealScalar;
70 typedef typename MatrixType::StorageIndex StorageIndex;
71 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
72 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
73 typedef Matrix<StorageIndex, 1,
74 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime),
RowMajor, 1,
75 EIGEN_SIZE_MIN_PREFER_FIXED(MaxColsAtCompileTime,MaxRowsAtCompileTime)> IntDiagSizeVectorType;
76 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationType;
77 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
78 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
79 typedef typename MatrixType::PlainObject PlainObject;
89 m_rows_transpositions(),
90 m_cols_transpositions(),
93 m_isInitialized(false),
94 m_usePrescribedThreshold(false) {}
104 m_hCoeffs((
std::min)(rows,cols)),
105 m_rows_transpositions((
std::min)(rows,cols)),
106 m_cols_transpositions((
std::min)(rows,cols)),
107 m_cols_permutation(cols),
109 m_isInitialized(false),
110 m_usePrescribedThreshold(false) {}
124 template<
typename InputType>
126 : m_qr(matrix.rows(), matrix.cols()),
127 m_hCoeffs((
std::min)(matrix.rows(), matrix.cols())),
128 m_rows_transpositions((
std::min)(matrix.rows(), matrix.cols())),
129 m_cols_transpositions((
std::min)(matrix.rows(), matrix.cols())),
130 m_cols_permutation(matrix.cols()),
131 m_temp(matrix.cols()),
132 m_isInitialized(false),
133 m_usePrescribedThreshold(false)
156 template<
typename Rhs>
160 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
166 MatrixQReturnType matrixQ(
void)
const;
172 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
176 template<
typename InputType>
182 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
183 return m_cols_permutation;
189 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
190 return m_rows_transpositions;
206 typename MatrixType::RealScalar absDeterminant()
const;
220 typename MatrixType::RealScalar logAbsDeterminant()
const;
231 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
232 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
234 for(Index i = 0; i < m_nonzero_pivots; ++i)
235 result += (abs(m_qr.coeff(i,i)) > premultiplied_threshold);
247 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
248 return cols() - rank();
260 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
261 return rank() == cols();
273 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
274 return rank() == rows();
285 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
286 return isInjective() && isSurjective();
296 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
300 inline Index rows()
const {
return m_qr.rows(); }
301 inline Index cols()
const {
return m_qr.cols(); }
307 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
328 m_usePrescribedThreshold =
true;
329 m_prescribedThreshold = threshold;
343 m_usePrescribedThreshold =
false;
353 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
354 return m_usePrescribedThreshold ? m_prescribedThreshold
369 eigen_assert(m_isInitialized &&
"LU is not initialized.");
370 return m_nonzero_pivots;
378 #ifndef EIGEN_PARSED_BY_DOXYGEN 379 template<
typename RhsType,
typename DstType>
381 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
386 static void check_template_parameters()
388 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
391 void computeInPlace();
394 HCoeffsType m_hCoeffs;
395 IntDiagSizeVectorType m_rows_transpositions;
396 IntDiagSizeVectorType m_cols_transpositions;
397 PermutationType m_cols_permutation;
398 RowVectorType m_temp;
399 bool m_isInitialized, m_usePrescribedThreshold;
400 RealScalar m_prescribedThreshold, m_maxpivot;
401 Index m_nonzero_pivots;
402 RealScalar m_precision;
406 template<
typename MatrixType>
410 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
411 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
412 return abs(m_qr.diagonal().prod());
415 template<
typename MatrixType>
418 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
419 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
420 return m_qr.diagonal().cwiseAbs().array().log().sum();
429 template<
typename MatrixType>
430 template<
typename InputType>
433 check_template_parameters();
442 template<
typename MatrixType>
446 Index rows = m_qr.rows();
447 Index cols = m_qr.cols();
448 Index size = (std::min)(rows,cols);
451 m_hCoeffs.resize(size);
457 m_rows_transpositions.resize(size);
458 m_cols_transpositions.resize(size);
459 Index number_of_transpositions = 0;
461 RealScalar biggest(0);
463 m_nonzero_pivots = size;
464 m_maxpivot = RealScalar(0);
466 for (Index k = 0; k < size; ++k)
468 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
469 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
470 typedef typename Scoring::result_type Score;
472 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
473 .unaryExpr(Scoring())
474 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
475 row_of_biggest_in_corner += k;
476 col_of_biggest_in_corner += k;
477 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
478 if(k==0) biggest = biggest_in_corner;
481 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
483 m_nonzero_pivots = k;
484 for(Index i = k; i < size; i++)
486 m_rows_transpositions.coeffRef(i) = i;
487 m_cols_transpositions.coeffRef(i) = i;
488 m_hCoeffs.coeffRef(i) = Scalar(0);
493 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
494 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
495 if(k != row_of_biggest_in_corner) {
496 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
497 ++number_of_transpositions;
499 if(k != col_of_biggest_in_corner) {
500 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
501 ++number_of_transpositions;
505 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
506 m_qr.coeffRef(k,k) = beta;
509 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
511 m_qr.bottomRightCorner(rows-k, cols-k-1)
512 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
515 m_cols_permutation.setIdentity(cols);
516 for(Index k = 0; k < size; ++k)
517 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
519 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
520 m_isInitialized =
true;
523 #ifndef EIGEN_PARSED_BY_DOXYGEN 524 template<
typename _MatrixType>
525 template<
typename RhsType,
typename DstType>
528 eigen_assert(rhs.rows() == rows());
529 const Index l_rank = rank();
539 typename RhsType::PlainObject c(rhs);
542 for (Index k = 0; k < l_rank; ++k)
544 Index remainingSize = rows()-k;
545 c.
row(k).swap(c.row(m_rows_transpositions.coeff(k)));
546 c.bottomRightCorner(remainingSize, rhs.cols())
547 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
548 m_hCoeffs.coeff(k), &temp.coeffRef(0));
551 m_qr.topLeftCorner(l_rank, l_rank)
552 .template triangularView<Upper>()
553 .solveInPlace(c.topRows(l_rank));
555 for(Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
556 for(Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
562 template<
typename DstXprType,
typename MatrixType,
typename Scalar>
563 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
567 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
569 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
579 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
580 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
584 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
585 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
586 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
588 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
589 const HCoeffsType& hCoeffs,
590 const IntDiagSizeVectorType& rowsTranspositions)
593 m_rowsTranspositions(rowsTranspositions)
596 template <
typename ResultType>
597 void evalTo(ResultType& result)
const 599 const Index rows = m_qr.rows();
600 WorkVectorType workspace(rows);
601 evalTo(result, workspace);
604 template <
typename ResultType>
605 void evalTo(ResultType& result, WorkVectorType& workspace)
const 611 const Index rows = m_qr.rows();
612 const Index cols = m_qr.cols();
613 const Index size = (std::min)(rows, cols);
614 workspace.resize(rows);
615 result.setIdentity(rows, rows);
616 for (Index k = size-1; k >= 0; k--)
618 result.block(k, k, rows-k, rows-k)
619 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
620 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
624 Index rows()
const {
return m_qr.rows(); }
625 Index cols()
const {
return m_qr.rows(); }
628 typename MatrixType::Nested m_qr;
629 typename HCoeffsType::Nested m_hCoeffs;
630 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
640 template<
typename MatrixType>
643 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
644 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
652 template<
typename Derived>
662 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:245
bool isInvertible() const
Definition: FullPivHouseholderQR.h:283
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:341
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: ForwardDeclarations.h:254
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:641
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:376
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:187
bool isInjective() const
Definition: FullPivHouseholderQR.h:258
Definition: StdDeque.h:58
bool isSurjective() const
Definition: FullPivHouseholderQR.h:271
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
RowXpr row(Index i)
Definition: DenseBase.h:797
Derived & derived()
Definition: EigenBase.h:44
Definition: EigenBase.h:28
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:102
Expression of the inverse of another expression.
Definition: Inverse.h:43
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivHouseholderQR.h:158
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:307
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:86
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:180
Definition: Constants.h:322
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:326
Definition: Eigen_Colamd.h:54
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:125
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:416
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:351
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:407
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:170
Index rank() const
Definition: FullPivHouseholderQR.h:228
Pseudo expression representing a solving operation.
Definition: Solve.h:62
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:654
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:294
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:367