11 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
18 template<
typename _MatrixType>
struct traits<FullPivHouseholderQR<_MatrixType> >
21 typedef MatrixXpr XprKind;
22 typedef SolverStorage StorageKind;
23 typedef int StorageIndex;
27 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType;
29 template<
typename MatrixType>
30 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
32 typedef typename MatrixType::PlainObject ReturnType;
61 :
public SolverBase<FullPivHouseholderQR<_MatrixType> >
65 typedef _MatrixType MatrixType;
71 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
72 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
74 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
75 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
76 typedef Matrix<StorageIndex, 1,
77 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime),
RowMajor, 1,
80 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
81 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
82 typedef typename MatrixType::PlainObject PlainObject;
92 m_rows_transpositions(),
93 m_cols_transpositions(),
96 m_isInitialized(false),
97 m_usePrescribedThreshold(false) {}
107 m_hCoeffs((std::min)(rows,cols)),
108 m_rows_transpositions((std::min)(rows,cols)),
109 m_cols_transpositions((std::min)(rows,cols)),
110 m_cols_permutation(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_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
132 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
133 m_cols_permutation(matrix.cols()),
134 m_temp(matrix.cols()),
135 m_isInitialized(false),
136 m_usePrescribedThreshold(false)
147 template<
typename InputType>
150 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
151 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
152 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
153 m_cols_permutation(matrix.cols()),
154 m_temp(matrix.cols()),
155 m_isInitialized(false),
156 m_usePrescribedThreshold(false)
161 #ifdef EIGEN_PARSED_BY_DOXYGEN
177 template<
typename Rhs>
184 MatrixQReturnType
matrixQ(
void)
const;
190 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
194 template<
typename InputType>
200 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
201 return m_cols_permutation;
207 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
208 return m_rows_transpositions;
249 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
250 RealScalar premultiplied_threshold =
abs(m_maxpivot) *
threshold();
252 for(
Index i = 0; i < m_nonzero_pivots; ++i)
253 result += (
abs(m_qr.coeff(i,i)) > premultiplied_threshold);
265 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
266 return cols() -
rank();
278 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
279 return rank() == cols();
291 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
292 return rank() == rows();
303 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
314 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
318 inline Index rows()
const {
return m_qr.rows(); }
319 inline Index cols()
const {
return m_qr.cols(); }
325 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
346 m_usePrescribedThreshold =
true;
361 m_usePrescribedThreshold =
false;
371 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
372 return m_usePrescribedThreshold ? m_prescribedThreshold
387 eigen_assert(m_isInitialized &&
"LU is not initialized.");
388 return m_nonzero_pivots;
396 #ifndef EIGEN_PARSED_BY_DOXYGEN
397 template<
typename RhsType,
typename DstType>
398 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
400 template<
bool Conjugate,
typename RhsType,
typename DstType>
401 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
406 static void check_template_parameters()
408 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
411 void computeInPlace();
414 HCoeffsType m_hCoeffs;
415 IntDiagSizeVectorType m_rows_transpositions;
416 IntDiagSizeVectorType m_cols_transpositions;
417 PermutationType m_cols_permutation;
418 RowVectorType m_temp;
419 bool m_isInitialized, m_usePrescribedThreshold;
420 RealScalar m_prescribedThreshold, m_maxpivot;
421 Index m_nonzero_pivots;
422 RealScalar m_precision;
426 template<
typename MatrixType>
430 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
431 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
432 return abs(m_qr.diagonal().prod());
435 template<
typename MatrixType>
438 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
439 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
440 return m_qr.diagonal().cwiseAbs().array().log().sum();
449 template<
typename MatrixType>
450 template<
typename InputType>
458 template<
typename MatrixType>
461 check_template_parameters();
464 Index rows = m_qr.rows();
465 Index cols = m_qr.cols();
466 Index size = (std::min)(rows,cols);
469 m_hCoeffs.resize(size);
475 m_rows_transpositions.resize(size);
476 m_cols_transpositions.resize(size);
477 Index number_of_transpositions = 0;
479 RealScalar biggest(0);
481 m_nonzero_pivots = size;
482 m_maxpivot = RealScalar(0);
484 for (
Index k = 0; k < size; ++k)
486 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
487 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
488 typedef typename Scoring::result_type Score;
490 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
491 .unaryExpr(Scoring())
492 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
493 row_of_biggest_in_corner += k;
494 col_of_biggest_in_corner += k;
495 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
496 if(k==0) biggest = biggest_in_corner;
499 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
501 m_nonzero_pivots = k;
502 for(
Index i = k; i < size; i++)
504 m_rows_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
505 m_cols_transpositions.coeffRef(i) = internal::convert_index<StorageIndex>(i);
506 m_hCoeffs.coeffRef(i) = Scalar(0);
511 m_rows_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(row_of_biggest_in_corner);
512 m_cols_transpositions.coeffRef(k) = internal::convert_index<StorageIndex>(col_of_biggest_in_corner);
513 if(k != row_of_biggest_in_corner) {
514 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
515 ++number_of_transpositions;
517 if(k != col_of_biggest_in_corner) {
518 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
519 ++number_of_transpositions;
523 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
524 m_qr.coeffRef(k,k) = beta;
527 if(
abs(beta) > m_maxpivot) m_maxpivot =
abs(beta);
529 m_qr.bottomRightCorner(rows-k, cols-k-1)
530 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
533 m_cols_permutation.setIdentity(cols);
534 for(
Index k = 0; k < size; ++k)
535 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
537 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
538 m_isInitialized =
true;
541 #ifndef EIGEN_PARSED_BY_DOXYGEN
542 template<
typename _MatrixType>
543 template<
typename RhsType,
typename DstType>
544 void FullPivHouseholderQR<_MatrixType>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
546 const Index l_rank = rank();
556 typename RhsType::PlainObject c(rhs);
558 Matrix<typename RhsType::Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
559 for (
Index k = 0; k < l_rank; ++k)
561 Index remainingSize = rows()-k;
562 c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
563 c.bottomRightCorner(remainingSize, rhs.cols())
564 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
565 m_hCoeffs.coeff(k), &temp.coeffRef(0));
568 m_qr.topLeftCorner(l_rank, l_rank)
569 .template triangularView<Upper>()
570 .solveInPlace(c.topRows(l_rank));
572 for(
Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
573 for(
Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
576 template<
typename _MatrixType>
577 template<
bool Conjugate,
typename RhsType,
typename DstType>
578 void FullPivHouseholderQR<_MatrixType>::_solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const
580 const Index l_rank = rank();
588 typename RhsType::PlainObject c(m_cols_permutation.transpose()*rhs);
590 m_qr.topLeftCorner(l_rank, l_rank)
591 .template triangularView<Upper>()
592 .transpose().template conjugateIf<Conjugate>()
593 .solveInPlace(c.topRows(l_rank));
595 dst.topRows(l_rank) = c.topRows(l_rank);
596 dst.bottomRows(rows()-l_rank).setZero();
598 Matrix<Scalar, 1, DstType::ColsAtCompileTime> temp(dst.cols());
599 const Index size = (std::min)(rows(), cols());
600 for (
Index k = size-1; k >= 0; --k)
602 Index remainingSize = rows()-k;
604 dst.bottomRightCorner(remainingSize, dst.cols())
605 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1).template conjugateIf<!Conjugate>(),
606 m_hCoeffs.template conjugateIf<Conjugate>().coeff(k), &temp.coeffRef(0));
608 dst.row(k).swap(dst.row(m_rows_transpositions.coeff(k)));
615 template<
typename DstXprType,
typename MatrixType>
616 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
618 typedef FullPivHouseholderQR<MatrixType> QrType;
619 typedef Inverse<QrType> SrcXprType;
620 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
622 dst = src.nestedExpression().
solve(MatrixType::Identity(src.rows(), src.cols()));
632 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
633 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
636 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
637 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
638 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
639 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
641 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
642 const HCoeffsType& hCoeffs,
643 const IntDiagSizeVectorType& rowsTranspositions)
646 m_rowsTranspositions(rowsTranspositions)
649 template <
typename ResultType>
650 void evalTo(ResultType& result)
const
652 const Index rows = m_qr.rows();
653 WorkVectorType workspace(rows);
654 evalTo(result, workspace);
657 template <
typename ResultType>
658 void evalTo(ResultType& result, WorkVectorType& workspace)
const
664 const Index rows = m_qr.rows();
665 const Index cols = m_qr.cols();
666 const Index size = (std::min)(rows, cols);
667 workspace.resize(rows);
668 result.setIdentity(rows, rows);
669 for (
Index k = size-1; k >= 0; k--)
671 result.block(k, k, rows-k, rows-k)
672 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
673 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
677 Index rows()
const {
return m_qr.rows(); }
678 Index cols()
const {
return m_qr.rows(); }
681 typename MatrixType::Nested m_qr;
682 typename HCoeffsType::Nested m_hCoeffs;
683 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
693 template<
typename MatrixType>
696 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
697 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
704 template<
typename Derived>
Householder rank-revealing QR decomposition of a matrix with full pivoting.
Definition: FullPivHouseholderQR.h:62
MatrixType::RealScalar absDeterminant() const
Definition: FullPivHouseholderQR.h:427
const Inverse< FullPivHouseholderQR > inverse() const
Definition: FullPivHouseholderQR.h:312
Index dimensionOfKernel() const
Definition: FullPivHouseholderQR.h:263
bool isInjective() const
Definition: FullPivHouseholderQR.h:276
const Solve< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
RealScalar maxPivot() const
Definition: FullPivHouseholderQR.h:394
const HCoeffsType & hCoeffs() const
Definition: FullPivHouseholderQR.h:325
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
Definition: FullPivHouseholderQR.h:344
const MatrixType & matrixQR() const
Definition: FullPivHouseholderQR.h:188
bool isSurjective() const
Definition: FullPivHouseholderQR.h:289
FullPivHouseholderQR & setThreshold(Default_t)
Definition: FullPivHouseholderQR.h:359
MatrixType::RealScalar logAbsDeterminant() const
Definition: FullPivHouseholderQR.h:436
const IntDiagSizeVectorType & rowsTranspositions() const
Definition: FullPivHouseholderQR.h:205
const PermutationType & colsPermutation() const
Definition: FullPivHouseholderQR.h:198
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
Definition: FullPivHouseholderQR.h:105
FullPivHouseholderQR(EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:148
MatrixQReturnType matrixQ(void) const
Definition: FullPivHouseholderQR.h:694
Index rank() const
Definition: FullPivHouseholderQR.h:246
FullPivHouseholderQR()
Default Constructor.
Definition: FullPivHouseholderQR.h:89
bool isInvertible() const
Definition: FullPivHouseholderQR.h:301
FullPivHouseholderQR(const EigenBase< InputType > &matrix)
Constructs a QR factorization from a given matrix.
Definition: FullPivHouseholderQR.h:128
Index nonzeroPivots() const
Definition: FullPivHouseholderQR.h:385
RealScalar threshold() const
Definition: FullPivHouseholderQR.h:369
Expression of the inverse of another expression.
Definition: Inverse.h:44
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Definition: FullPivHouseholderQR.h:706
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
Pseudo expression representing a solving operation.
Definition: Solve.h:63
A base class for matrix decomposition and solvers.
Definition: SolverBase.h:69
FullPivHouseholderQR< _MatrixType > & derived()
Definition: EigenBase.h:46
const Solve< Derived, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: SolverBase.h:106
@ RowMajor
Definition: Constants.h:321
Namespace containing all symbols from the Eigen library.
Definition: Core:134
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
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:30
Derived & derived()
Definition: EigenBase.h:46
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:39
Holds information about the various numeric (i.e. scalar) types allowed by Eigen.
Definition: NumTraits.h:213