11 #ifndef EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_COLPIVOTINGHOUSEHOLDERQR_H
14 #include "./InternalHeaderCheck.h"
19 template<
typename MatrixType_>
struct traits<ColPivHouseholderQR<MatrixType_> >
22 typedef MatrixXpr XprKind;
23 typedef SolverStorage StorageKind;
24 typedef int StorageIndex;
54 :
public SolverBase<ColPivHouseholderQR<MatrixType_> >
58 typedef MatrixType_ MatrixType;
64 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
65 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
67 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
69 typedef typename internal::plain_row_type<MatrixType, Index>::type IntRowVectorType;
70 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
71 typedef typename internal::plain_row_type<MatrixType, RealScalar>::type RealRowVectorType;
73 typedef typename MatrixType::PlainObject PlainObject;
77 typedef typename PermutationType::StorageIndex PermIndexType;
91 m_colsTranspositions(),
95 m_isInitialized(false),
96 m_usePrescribedThreshold(false) {}
106 m_hCoeffs((std::min)(rows,cols)),
107 m_colsPermutation(PermIndexType(cols)),
108 m_colsTranspositions(cols),
110 m_colNormsUpdated(cols),
111 m_colNormsDirect(cols),
112 m_isInitialized(false),
113 m_usePrescribedThreshold(false) {}
127 template<
typename InputType>
129 : m_qr(matrix.rows(), matrix.cols()),
130 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
131 m_colsPermutation(PermIndexType(matrix.cols())),
132 m_colsTranspositions(matrix.cols()),
133 m_temp(matrix.cols()),
134 m_colNormsUpdated(matrix.cols()),
135 m_colNormsDirect(matrix.cols()),
136 m_isInitialized(false),
137 m_usePrescribedThreshold(false)
148 template<
typename InputType>
151 m_hCoeffs((std::min)(matrix.rows(),matrix.cols())),
152 m_colsPermutation(PermIndexType(matrix.cols())),
153 m_colsTranspositions(matrix.cols()),
154 m_temp(matrix.cols()),
155 m_colNormsUpdated(matrix.cols()),
156 m_colNormsDirect(matrix.cols()),
157 m_isInitialized(false),
158 m_usePrescribedThreshold(false)
163 #ifdef EIGEN_PARSED_BY_DOXYGEN
178 template<
typename Rhs>
193 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
208 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
212 template<
typename InputType>
218 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
219 return m_colsPermutation;
260 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
261 RealScalar premultiplied_threshold =
abs(m_maxpivot) *
threshold();
263 for(
Index i = 0; i < m_nonzero_pivots; ++i)
264 result += (
abs(m_qr.coeff(i,i)) > premultiplied_threshold);
276 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
277 return cols() -
rank();
289 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
290 return rank() == cols();
302 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
303 return rank() == rows();
314 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
325 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
329 inline Index rows()
const {
return m_qr.rows(); }
330 inline Index cols()
const {
return m_qr.cols(); }
336 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
357 m_usePrescribedThreshold =
true;
372 m_usePrescribedThreshold =
false;
382 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
383 return m_usePrescribedThreshold ? m_prescribedThreshold
398 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
399 return m_nonzero_pivots;
415 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
419 #ifndef EIGEN_PARSED_BY_DOXYGEN
420 template<
typename RhsType,
typename DstType>
421 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
423 template<
bool Conjugate,
typename RhsType,
typename DstType>
424 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
431 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar)
433 void computeInPlace();
436 HCoeffsType m_hCoeffs;
437 PermutationType m_colsPermutation;
438 IntRowVectorType m_colsTranspositions;
439 RowVectorType m_temp;
440 RealRowVectorType m_colNormsUpdated;
441 RealRowVectorType m_colNormsDirect;
442 bool m_isInitialized, m_usePrescribedThreshold;
443 RealScalar m_prescribedThreshold, m_maxpivot;
444 Index m_nonzero_pivots;
448 template<
typename MatrixType>
452 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
453 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
454 return abs(m_qr.diagonal().prod());
457 template<
typename MatrixType>
460 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
461 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
462 return m_qr.diagonal().cwiseAbs().array().log().sum();
471 template<
typename MatrixType>
472 template<
typename InputType>
480 template<
typename MatrixType>
488 Index rows = m_qr.rows();
489 Index cols = m_qr.cols();
490 Index size = m_qr.diagonalSize();
492 m_hCoeffs.resize(size);
496 m_colsTranspositions.resize(m_qr.cols());
497 Index number_of_transpositions = 0;
499 m_colNormsUpdated.resize(cols);
500 m_colNormsDirect.resize(cols);
501 for (
Index k = 0; k < cols; ++k) {
504 m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
505 m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
508 RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
509 RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
511 m_nonzero_pivots = size;
512 m_maxpivot = RealScalar(0);
514 for(Index k = 0; k < size; ++k)
517 Index biggest_col_index;
518 RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
519 biggest_col_index += k;
523 if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
524 m_nonzero_pivots = k;
527 m_colsTranspositions.coeffRef(k) = biggest_col_index;
528 if(k != biggest_col_index) {
529 m_qr.col(k).swap(m_qr.col(biggest_col_index));
530 std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
531 std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
532 ++number_of_transpositions;
537 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
540 m_qr.coeffRef(k,k) = beta;
543 if(abs(beta) > m_maxpivot) m_maxpivot = abs(beta);
546 m_qr.bottomRightCorner(rows-k, cols-k-1)
547 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
550 for (Index j = k + 1; j < cols; ++j) {
555 if (!numext::is_exactly_zero(m_colNormsUpdated.coeffRef(j))) {
556 RealScalar temp = abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
557 temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
558 temp = temp < RealScalar(0) ? RealScalar(0) : temp;
559 RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
560 m_colNormsDirect.coeffRef(j));
561 if (temp2 <= norm_downdate_threshold) {
564 m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
565 m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
567 m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
573 m_colsPermutation.setIdentity(PermIndexType(cols));
574 for(PermIndexType k = 0; k < size; ++k)
575 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
577 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
578 m_isInitialized =
true;
581 #ifndef EIGEN_PARSED_BY_DOXYGEN
582 template<
typename MatrixType_>
583 template<
typename RhsType,
typename DstType>
584 void ColPivHouseholderQR<MatrixType_>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
586 const Index nonzero_pivots = nonzeroPivots();
588 if(nonzero_pivots == 0)
594 typename RhsType::PlainObject c(rhs);
596 c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
598 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
599 .template triangularView<Upper>()
600 .solveInPlace(c.topRows(nonzero_pivots));
602 for(Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
603 for(Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
606 template<
typename MatrixType_>
607 template<
bool Conjugate,
typename RhsType,
typename DstType>
608 void ColPivHouseholderQR<MatrixType_>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
610 const Index nonzero_pivots = nonzeroPivots();
612 if(nonzero_pivots == 0)
618 typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs);
620 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
621 .template triangularView<Upper>()
622 .transpose().template conjugateIf<Conjugate>()
623 .solveInPlace(c.topRows(nonzero_pivots));
625 dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
626 dst.bottomRows(rows()-nonzero_pivots).setZero();
628 dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).
template conjugateIf<!Conjugate>() );
634 template<
typename DstXprType,
typename MatrixType>
635 struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
637 typedef ColPivHouseholderQR<MatrixType> QrType;
638 typedef Inverse<QrType> SrcXprType;
639 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
641 dst = src.nestedExpression().
solve(MatrixType::Identity(src.rows(), src.cols()));
650 template<
typename MatrixType>
651 typename ColPivHouseholderQR<MatrixType>::HouseholderSequenceType ColPivHouseholderQR<MatrixType>
652 ::householderQ()
const
654 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
662 template<
typename Derived>
Householder rank-revealing QR decomposition of a matrix with column-pivoting.
Definition: ColPivHouseholderQR.h:55
const HCoeffsType & hCoeffs() const
Definition: ColPivHouseholderQR.h:336
ColPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:149
Index nonzeroPivots() const
Definition: ColPivHouseholderQR.h:396
bool isInvertible() const
Definition: ColPivHouseholderQR.h:312
ColPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: ColPivHouseholderQR.h:128
HouseholderSequenceType householderQ() const
Definition: ColPivHouseholderQR.h:652
const MatrixType & matrixQR() const
Definition: ColPivHouseholderQR.h:191
ComputationInfo info() const
Reports whether the QR factorization was successful.
Definition: ColPivHouseholderQR.h:413
Index rank() const
Definition: ColPivHouseholderQR.h:257
const Inverse< ColPivHouseholderQR > inverse() const
Definition: ColPivHouseholderQR.h:323
ColPivHouseholderQR()
Default Constructor.
Definition: ColPivHouseholderQR.h:87
ColPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: ColPivHouseholderQR.h:104
RealScalar threshold() const
Definition: ColPivHouseholderQR.h:380
const PermutationType & colsPermutation() const
Definition: ColPivHouseholderQR.h:216
RealScalar maxPivot() const
Definition: ColPivHouseholderQR.h:405
const MatrixType & matrixR() const
Definition: ColPivHouseholderQR.h:206
bool isSurjective() const
Definition: ColPivHouseholderQR.h:300
MatrixType::RealScalar absDeterminant() const
Definition: ColPivHouseholderQR.h:449
const Solve< ColPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
ColPivHouseholderQR & setThreshold(Default_t)
Definition: ColPivHouseholderQR.h:370
Index dimensionOfKernel() const
Definition: ColPivHouseholderQR.h:274
ColPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: ColPivHouseholderQR.h:355
bool isInjective() const
Definition: ColPivHouseholderQR.h:287
MatrixType::RealScalar logAbsDeterminant() const
Definition: ColPivHouseholderQR.h:458
Complete orthogonal decomposition (COD) of a matrix.
Definition: CompleteOrthogonalDecomposition.h:54
Sequence of Householder reflections acting on subspaces with decreasing size.
Definition: HouseholderSequence.h:123
Expression of the inverse of another expression.
Definition: Inverse.h:46
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:52
Pseudo expression representing a solving operation.
Definition: Solve.h:65
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:71
ColPivHouseholderQR< MatrixType_ > & derived()
Definition: EigenBase.h:48
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:107
ComputationInfo
Definition: Constants.h:442
@ Success
Definition: Constants.h:444
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_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