11 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
14 #include "./InternalHeaderCheck.h"
20 template<
typename MatrixType_>
struct traits<FullPivHouseholderQR<MatrixType_> >
23 typedef MatrixXpr XprKind;
24 typedef SolverStorage StorageKind;
25 typedef int StorageIndex;
29 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType;
31 template<
typename MatrixType>
32 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
34 typedef typename MatrixType::PlainObject ReturnType;
63 :
public SolverBase<FullPivHouseholderQR<MatrixType_> >
67 typedef MatrixType_ MatrixType;
73 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
74 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
76 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
77 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
78 typedef Matrix<StorageIndex, 1,
79 internal::min_size_prefer_dynamic(ColsAtCompileTime,RowsAtCompileTime),
RowMajor, 1,
82 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
83 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
84 typedef typename MatrixType::PlainObject PlainObject;
94 m_rows_transpositions(),
95 m_cols_transpositions(),
98 m_isInitialized(false),
99 m_usePrescribedThreshold(false) {}
109 m_hCoeffs((std::min)(rows,cols)),
110 m_rows_transpositions((std::min)(rows,cols)),
111 m_cols_transpositions((std::min)(rows,cols)),
112 m_cols_permutation(cols),
114 m_isInitialized(false),
115 m_usePrescribedThreshold(false) {}
129 template<
typename InputType>
131 : m_qr(matrix.rows(), matrix.cols()),
132 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
133 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
134 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
135 m_cols_permutation(matrix.cols()),
136 m_temp(matrix.cols()),
137 m_isInitialized(false),
138 m_usePrescribedThreshold(false)
149 template<
typename InputType>
152 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
153 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
154 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
155 m_cols_permutation(matrix.cols()),
156 m_temp(matrix.cols()),
157 m_isInitialized(false),
158 m_usePrescribedThreshold(false)
163 #ifdef EIGEN_PARSED_BY_DOXYGEN
179 template<
typename Rhs>
186 MatrixQReturnType
matrixQ(
void)
const;
192 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
196 template<
typename InputType>
202 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
203 return m_cols_permutation;
209 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
210 return m_rows_transpositions;
251 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
252 RealScalar premultiplied_threshold =
abs(m_maxpivot) *
threshold();
254 for(
Index i = 0; i < m_nonzero_pivots; ++i)
255 result += (
abs(m_qr.coeff(i,i)) > premultiplied_threshold);
267 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
268 return cols() -
rank();
280 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
281 return rank() == cols();
293 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
294 return rank() == rows();
305 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
316 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
320 inline Index rows()
const {
return m_qr.rows(); }
321 inline Index cols()
const {
return m_qr.cols(); }
327 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
348 m_usePrescribedThreshold =
true;
363 m_usePrescribedThreshold =
false;
373 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
374 return m_usePrescribedThreshold ? m_prescribedThreshold
389 eigen_assert(m_isInitialized &&
"LU is not initialized.");
390 return m_nonzero_pivots;
398 #ifndef EIGEN_PARSED_BY_DOXYGEN
399 template<
typename RhsType,
typename DstType>
400 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
402 template<
bool Conjugate,
typename RhsType,
typename DstType>
403 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
408 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
410 void computeInPlace();
413 HCoeffsType m_hCoeffs;
414 IntDiagSizeVectorType m_rows_transpositions;
415 IntDiagSizeVectorType m_cols_transpositions;
416 PermutationType m_cols_permutation;
417 RowVectorType m_temp;
418 bool m_isInitialized, m_usePrescribedThreshold;
419 RealScalar m_prescribedThreshold, m_maxpivot;
420 Index m_nonzero_pivots;
421 RealScalar m_precision;
425 template<
typename MatrixType>
429 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
430 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
431 return abs(m_qr.diagonal().prod());
434 template<
typename MatrixType>
437 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
438 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
439 return m_qr.diagonal().cwiseAbs().array().log().sum();
448 template<
typename MatrixType>
449 template<
typename InputType>
457 template<
typename MatrixType>
461 Index rows = m_qr.rows();
462 Index cols = m_qr.cols();
463 Index size = (std::min)(rows,cols);
466 m_hCoeffs.resize(size);
472 m_rows_transpositions.resize(size);
473 m_cols_transpositions.resize(size);
474 Index number_of_transpositions = 0;
476 RealScalar biggest(0);
478 m_nonzero_pivots = size;
479 m_maxpivot = RealScalar(0);
481 for (
Index k = 0; k < size; ++k)
483 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
484 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
485 typedef typename Scoring::result_type Score;
487 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
488 .unaryExpr(Scoring())
489 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
490 row_of_biggest_in_corner += k;
491 col_of_biggest_in_corner += k;
492 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
493 if(k==0) biggest = biggest_in_corner;
496 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
498 m_nonzero_pivots = k;
499 for(
Index i = k; i < size; i++)
501 m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
502 m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
503 m_hCoeffs.coeffRef(i) = Scalar(0);
508 m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner);
509 m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner);
510 if(k != row_of_biggest_in_corner) {
511 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
512 ++number_of_transpositions;
514 if(k != col_of_biggest_in_corner) {
515 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
516 ++number_of_transpositions;
520 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
521 m_qr.coeffRef(k,k) = beta;
524 if(
abs(beta) > m_maxpivot) m_maxpivot =
abs(beta);
526 m_qr.bottomRightCorner(rows-k, cols-k-1)
527 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
530 m_cols_permutation.setIdentity(cols);
531 for(
Index k = 0; k < size; ++k)
532 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
534 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
535 m_isInitialized =
true;
538 #ifndef EIGEN_PARSED_BY_DOXYGEN
539 template<
typename MatrixType_>
540 template<
typename RhsType,
typename DstType>
541 void FullPivHouseholderQR<MatrixType_>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
543 const Index l_rank = rank();
553 typename RhsType::PlainObject c(rhs);
555 Matrix<typename RhsType::Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
556 for (
Index k = 0; k < l_rank; ++k)
558 Index remainingSize = rows()-k;
559 c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
560 c.bottomRightCorner(remainingSize, rhs.cols())
561 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
562 m_hCoeffs.coeff(k), &temp.coeffRef(0));
565 m_qr.topLeftCorner(l_rank, l_rank)
566 .template triangularView<Upper>()
567 .solveInPlace(c.topRows(l_rank));
569 for(
Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
570 for(
Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
573 template<
typename MatrixType_>
574 template<
bool Conjugate,
typename RhsType,
typename DstType>
575 void FullPivHouseholderQR<MatrixType_>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
577 const Index l_rank = rank();
585 typename RhsType::PlainObject c(m_cols_permutation.transpose()*rhs);
587 m_qr.topLeftCorner(l_rank, l_rank)
588 .template triangularView<Upper>()
589 .transpose().template conjugateIf<Conjugate>()
590 .solveInPlace(c.topRows(l_rank));
592 dst.topRows(l_rank) = c.topRows(l_rank);
593 dst.bottomRows(rows()-l_rank).setZero();
595 Matrix<Scalar, 1, DstType::ColsAtCompileTime> temp(dst.cols());
596 const Index size = (std::min)(rows(), cols());
597 for (
Index k = size-1; k >= 0; --k)
599 Index remainingSize = rows()-k;
601 dst.bottomRightCorner(remainingSize, dst.cols())
602 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1).template conjugateIf<!Conjugate>(),
603 m_hCoeffs.template conjugateIf<Conjugate>().coeff(k), &temp.coeffRef(0));
605 dst.row(k).swap(dst.row(m_rows_transpositions.coeff(k)));
612 template<
typename DstXprType,
typename MatrixType>
613 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
615 typedef FullPivHouseholderQR<MatrixType> QrType;
616 typedef Inverse<QrType> SrcXprType;
617 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
619 dst = src.nestedExpression().
solve(MatrixType::Identity(src.rows(), src.cols()));
629 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
630 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
633 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
634 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
635 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
636 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
638 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
639 const HCoeffsType& hCoeffs,
640 const IntDiagSizeVectorType& rowsTranspositions)
643 m_rowsTranspositions(rowsTranspositions)
646 template <
typename ResultType>
647 void evalTo(ResultType& result)
const
649 const Index rows = m_qr.rows();
650 WorkVectorType workspace(rows);
651 evalTo(result, workspace);
654 template <
typename ResultType>
655 void evalTo(ResultType& result, WorkVectorType& workspace)
const
661 const Index rows = m_qr.rows();
662 const Index cols = m_qr.cols();
663 const Index size = (std::min)(rows, cols);
664 workspace.resize(rows);
665 result.setIdentity(rows, rows);
666 for (
Index k = size-1; k >= 0; k--)
668 result.block(k, k, rows-k, rows-k)
669 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
670 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
674 Index rows()
const {
return m_qr.rows(); }
675 Index cols()
const {
return m_qr.rows(); }
678 typename MatrixType::Nested m_qr;
679 typename HCoeffsType::Nested m_hCoeffs;
680 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
690 template<
typename MatrixType>
693 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
694 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
701 template<
typename Derived>
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: FullPivHouseholderQR.h:64
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:426
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:130
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:387
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:265
bool isInvertible() const
Definition: FullPivHouseholderQR.h:303
FullPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:150
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:107
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:91
Index rank() const
Definition: FullPivHouseholderQR.h:248
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:327
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:200
bool isInjective() const
Definition: FullPivHouseholderQR.h:278
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:435
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:190
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:396
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:361
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:346
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:207
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:691
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:371
bool isSurjective() const
Definition: FullPivHouseholderQR.h:291
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:314
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:703
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:182
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
FullPivHouseholderQR< MatrixType_ > & derived()
Definition: EigenBase.h:48
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:107
@ RowMajor
Definition: Constants.h:323
Namespace containing all symbols from the Eigen library.
Definition: Core:139
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:59
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
Definition: EigenBase.h:32
Derived & derived()
Definition: EigenBase.h:48
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:41
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:231