370 m_usePrescribedThreshold =
false;
380 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
381 return m_usePrescribedThreshold ? m_prescribedThreshold
384 : NumTraits<Scalar>::epsilon() * RealScalar(m_qr.diagonalSize());
396 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
397 return m_nonzero_pivots;
413 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
417 #ifndef EIGEN_PARSED_BY_DOXYGEN
418 template<
typename RhsType,
typename DstType>
419 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
421 template<
bool Conjugate,
typename RhsType,
typename DstType>
422 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
429 static void check_template_parameters()
431 EIGEN_STATIC_ASSERT_NON_INTEGER(
Scalar);
434 void computeInPlace();
437 HCoeffsType m_hCoeffs;
438 PermutationType m_colsPermutation;
439 IntRowVectorType m_colsTranspositions;
440 RowVectorType m_temp;
441 RealRowVectorType m_colNormsUpdated;
442 RealRowVectorType m_colNormsDirect;
443 bool m_isInitialized, m_usePrescribedThreshold;
444 RealScalar m_prescribedThreshold, m_maxpivot;
445 Index m_nonzero_pivots;
449template<
typename MatrixType>
453 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
454 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
455 return abs(m_qr.diagonal().prod());
458template<
typename MatrixType>
461 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
462 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
463 return m_qr.diagonal().cwiseAbs().array().log().sum();
472template<
typename MatrixType>
473template<
typename InputType>
481template<
typename MatrixType>
482void ColPivHouseholderQR<MatrixType>::computeInPlace()
484 check_template_parameters();
487 eigen_assert(m_qr.cols()<=NumTraits<int>::highest());
491 Index rows = m_qr.rows();
492 Index cols = m_qr.cols();
493 Index size = m_qr.diagonalSize();
495 m_hCoeffs.resize(size);
499 m_colsTranspositions.resize(m_qr.cols());
500 Index number_of_transpositions = 0;
502 m_colNormsUpdated.resize(cols);
503 m_colNormsDirect.resize(cols);
504 for (
Index k = 0; k < cols; ++k) {
507 m_colNormsDirect.coeffRef(k) = m_qr.col(k).norm();
508 m_colNormsUpdated.coeffRef(k) = m_colNormsDirect.coeffRef(k);
511 RealScalar threshold_helper = numext::abs2<RealScalar>(m_colNormsUpdated.maxCoeff() * NumTraits<RealScalar>::epsilon()) / RealScalar(rows);
512 RealScalar norm_downdate_threshold = numext::sqrt(NumTraits<RealScalar>::epsilon());
514 m_nonzero_pivots = size;
515 m_maxpivot = RealScalar(0);
517 for(Index k = 0; k < size; ++k)
520 Index biggest_col_index;
521 RealScalar biggest_col_sq_norm = numext::abs2(m_colNormsUpdated.tail(cols-k).maxCoeff(&biggest_col_index));
522 biggest_col_index += k;
526 if(m_nonzero_pivots==size && biggest_col_sq_norm < threshold_helper * RealScalar(rows-k))
527 m_nonzero_pivots = k;
530 m_colsTranspositions.coeffRef(k) = biggest_col_index;
531 if(k != biggest_col_index) {
532 m_qr.col(k).swap(m_qr.col(biggest_col_index));
533 std::swap(m_colNormsUpdated.coeffRef(k), m_colNormsUpdated.coeffRef(biggest_col_index));
534 std::swap(m_colNormsDirect.coeffRef(k), m_colNormsDirect.coeffRef(biggest_col_index));
535 ++number_of_transpositions;
540 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
543 m_qr.coeffRef(k,k) = beta;
546 if(
abs(beta) > m_maxpivot) m_maxpivot =
abs(beta);
549 m_qr.bottomRightCorner(rows-k, cols-k-1)
550 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
553 for (
Index j = k + 1; j < cols; ++j) {
558 if (m_colNormsUpdated.coeffRef(j) != RealScalar(0)) {
559 RealScalar temp =
abs(m_qr.coeffRef(k, j)) / m_colNormsUpdated.coeffRef(j);
560 temp = (RealScalar(1) + temp) * (RealScalar(1) - temp);
561 temp = temp < RealScalar(0) ? RealScalar(0) : temp;
562 RealScalar temp2 = temp * numext::abs2<RealScalar>(m_colNormsUpdated.coeffRef(j) /
563 m_colNormsDirect.coeffRef(j));
564 if (temp2 <= norm_downdate_threshold) {
567 m_colNormsDirect.coeffRef(j) = m_qr.col(j).tail(rows - k - 1).norm();
568 m_colNormsUpdated.coeffRef(j) = m_colNormsDirect.coeffRef(j);
570 m_colNormsUpdated.coeffRef(j) *= numext::sqrt(temp);
576 m_colsPermutation.setIdentity(PermIndexType(cols));
577 for(PermIndexType k = 0; k < size; ++k)
578 m_colsPermutation.applyTranspositionOnTheRight(k, PermIndexType(m_colsTranspositions.coeff(k)));
580 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
581 m_isInitialized =
true;
584#ifndef EIGEN_PARSED_BY_DOXYGEN
585template<
typename _MatrixType>
586template<
typename RhsType,
typename DstType>
589 const Index nonzero_pivots = nonzeroPivots();
591 if(nonzero_pivots == 0)
597 typename RhsType::PlainObject c(rhs);
599 c.applyOnTheLeft(householderQ().setLength(nonzero_pivots).adjoint() );
601 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
602 .template triangularView<Upper>()
603 .solveInPlace(c.topRows(nonzero_pivots));
605 for(
Index i = 0; i < nonzero_pivots; ++i) dst.row(m_colsPermutation.indices().coeff(i)) = c.row(i);
606 for(
Index i = nonzero_pivots; i < cols(); ++i) dst.row(m_colsPermutation.indices().coeff(i)).setZero();
609template<
typename _MatrixType>
610template<
bool Conjugate,
typename RhsType,
typename DstType>
613 const Index nonzero_pivots = nonzeroPivots();
615 if(nonzero_pivots == 0)
621 typename RhsType::PlainObject c(m_colsPermutation.transpose()*rhs);
623 m_qr.topLeftCorner(nonzero_pivots, nonzero_pivots)
624 .template triangularView<Upper>()
625 .transpose().template conjugateIf<Conjugate>()
626 .solveInPlace(c.topRows(nonzero_pivots));
628 dst.topRows(nonzero_pivots) = c.topRows(nonzero_pivots);
629 dst.bottomRows(rows()-nonzero_pivots).setZero();
631 dst.applyOnTheLeft(householderQ().setLength(nonzero_pivots).
template conjugateIf<!Conjugate>() );
637template<
typename DstXprType,
typename MatrixType>
638struct Assignment<DstXprType, Inverse<ColPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename ColPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
640 typedef ColPivHouseholderQR<MatrixType> QrType;
641 typedef Inverse<QrType> SrcXprType;
642 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
644 dst = src.nestedExpression().
solve(MatrixType::Identity(src.rows(), src.cols()));
653template<
typename MatrixType>
657 eigen_assert(m_isInitialized &&
"ColPivHouseholderQR is not initialized.");
658 return HouseholderSequenceType(m_qr, m_hCoeffs.conjugate());
665template<
typename Derived>