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;
61 typedef _MatrixType MatrixType;
63 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
64 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
65 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
66 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
68 typedef typename MatrixType::Scalar Scalar;
69 typedef typename MatrixType::RealScalar RealScalar;
71 typedef typename MatrixType::StorageIndex StorageIndex;
72 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
73 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
74 typedef Matrix<StorageIndex, 1,
75 EIGEN_SIZE_MIN_PREFER_DYNAMIC(ColsAtCompileTime,RowsAtCompileTime),
RowMajor, 1,
78 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
79 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
80 typedef typename MatrixType::PlainObject PlainObject;
90 m_rows_transpositions(),
91 m_cols_transpositions(),
94 m_isInitialized(false),
95 m_usePrescribedThreshold(false) {}
105 m_hCoeffs((std::min)(rows,cols)),
106 m_rows_transpositions((std::min)(rows,cols)),
107 m_cols_transpositions((std::min)(rows,cols)),
108 m_cols_permutation(cols),
110 m_isInitialized(false),
111 m_usePrescribedThreshold(false) {}
125 template<
typename InputType>
127 : m_qr(matrix.rows(), matrix.cols()),
128 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
129 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
130 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
131 m_cols_permutation(matrix.cols()),
132 m_temp(matrix.cols()),
133 m_isInitialized(false),
134 m_usePrescribedThreshold(false)
145 template<
typename InputType>
147 : m_qr(matrix.derived()),
148 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
149 m_rows_transpositions((std::min)(matrix.rows(), matrix.cols())),
150 m_cols_transpositions((std::min)(matrix.rows(), matrix.cols())),
151 m_cols_permutation(matrix.cols()),
152 m_temp(matrix.cols()),
153 m_isInitialized(false),
154 m_usePrescribedThreshold(false)
174 template<
typename Rhs>
178 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
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>
399 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
404 static void check_template_parameters()
406 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
409 void computeInPlace();
412 HCoeffsType m_hCoeffs;
413 IntDiagSizeVectorType m_rows_transpositions;
414 IntDiagSizeVectorType m_cols_transpositions;
415 PermutationType m_cols_permutation;
416 RowVectorType m_temp;
417 bool m_isInitialized, m_usePrescribedThreshold;
418 RealScalar m_prescribedThreshold, m_maxpivot;
419 Index m_nonzero_pivots;
420 RealScalar m_precision;
424 template<
typename MatrixType>
428 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
429 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
430 return abs(m_qr.diagonal().prod());
433 template<
typename MatrixType>
436 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
437 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
438 return m_qr.diagonal().cwiseAbs().array().log().sum();
447 template<
typename MatrixType>
448 template<
typename InputType>
456 template<
typename MatrixType>
459 check_template_parameters();
462 Index rows = m_qr.rows();
463 Index cols = m_qr.cols();
464 Index size = (std::min)(rows,cols);
467 m_hCoeffs.resize(size);
473 m_rows_transpositions.resize(size);
474 m_cols_transpositions.resize(size);
475 Index number_of_transpositions = 0;
477 RealScalar biggest(0);
479 m_nonzero_pivots = size;
480 m_maxpivot = RealScalar(0);
482 for (
Index k = 0; k < size; ++k)
484 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
485 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
486 typedef typename Scoring::result_type Score;
488 Score score = m_qr.bottomRightCorner(rows-k, cols-k)
489 .unaryExpr(Scoring())
490 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
491 row_of_biggest_in_corner += k;
492 col_of_biggest_in_corner += k;
493 RealScalar biggest_in_corner = internal::abs_knowing_score<Scalar>()(m_qr(row_of_biggest_in_corner, col_of_biggest_in_corner), score);
494 if(k==0) biggest = biggest_in_corner;
497 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
499 m_nonzero_pivots = k;
500 for(
Index i = k; i < size; i++)
502 m_rows_transpositions.coeffRef(i) = i;
503 m_cols_transpositions.coeffRef(i) = i;
504 m_hCoeffs.coeffRef(i) = Scalar(0);
509 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
510 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
511 if(k != row_of_biggest_in_corner) {
512 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
513 ++number_of_transpositions;
515 if(k != col_of_biggest_in_corner) {
516 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
517 ++number_of_transpositions;
521 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
522 m_qr.coeffRef(k,k) = beta;
525 if(
abs(beta) > m_maxpivot) m_maxpivot =
abs(beta);
527 m_qr.bottomRightCorner(rows-k, cols-k-1)
528 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
531 m_cols_permutation.setIdentity(cols);
532 for(
Index k = 0; k < size; ++k)
533 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
535 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
536 m_isInitialized =
true;
539 #ifndef EIGEN_PARSED_BY_DOXYGEN
540 template<
typename _MatrixType>
541 template<
typename RhsType,
typename DstType>
542 void FullPivHouseholderQR<_MatrixType>::_solve_impl(
const RhsType &rhs, DstType &dst)
const
544 eigen_assert(rhs.rows() == rows());
545 const Index l_rank = rank();
555 typename RhsType::PlainObject c(rhs);
557 Matrix<Scalar,1,RhsType::ColsAtCompileTime> temp(rhs.cols());
558 for (
Index k = 0; k < l_rank; ++k)
560 Index remainingSize = rows()-k;
561 c.row(k).swap(c.row(m_rows_transpositions.coeff(k)));
562 c.bottomRightCorner(remainingSize, rhs.cols())
563 .applyHouseholderOnTheLeft(m_qr.col(k).tail(remainingSize-1),
564 m_hCoeffs.coeff(k), &temp.coeffRef(0));
567 m_qr.topLeftCorner(l_rank, l_rank)
568 .template triangularView<Upper>()
569 .solveInPlace(c.topRows(l_rank));
571 for(
Index i = 0; i < l_rank; ++i) dst.row(m_cols_permutation.indices().coeff(i)) = c.row(i);
572 for(
Index i = l_rank; i < cols(); ++i) dst.row(m_cols_permutation.indices().coeff(i)).setZero();
578 template<
typename DstXprType,
typename MatrixType>
579 struct Assignment<DstXprType, Inverse<FullPivHouseholderQR<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename FullPivHouseholderQR<MatrixType>::Scalar>, Dense2Dense>
581 typedef FullPivHouseholderQR<MatrixType> QrType;
582 typedef Inverse<QrType> SrcXprType;
583 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename QrType::Scalar> &)
585 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
595 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
596 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
599 typedef typename FullPivHouseholderQR<MatrixType>::IntDiagSizeVectorType IntDiagSizeVectorType;
600 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
601 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
602 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
604 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
605 const HCoeffsType& hCoeffs,
606 const IntDiagSizeVectorType& rowsTranspositions)
609 m_rowsTranspositions(rowsTranspositions)
612 template <
typename ResultType>
613 void evalTo(ResultType& result)
const
615 const Index rows = m_qr.rows();
616 WorkVectorType workspace(rows);
617 evalTo(result, workspace);
620 template <
typename ResultType>
621 void evalTo(ResultType& result, WorkVectorType& workspace)
const
627 const Index rows = m_qr.rows();
628 const Index cols = m_qr.cols();
629 const Index size = (std::min)(rows, cols);
630 workspace.resize(rows);
631 result.setIdentity(rows, rows);
632 for (
Index k = size-1; k >= 0; k--)
634 result.block(k, k, rows-k, rows-k)
635 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
636 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
640 Index rows()
const {
return m_qr.rows(); }
641 Index cols()
const {
return m_qr.rows(); }
644 typename MatrixType::Nested m_qr;
645 typename HCoeffsType::Nested m_hCoeffs;
646 typename IntDiagSizeVectorType::Nested m_rowsTranspositions;
656 template<
typename MatrixType>
659 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
660 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
667 template<
typename Derived>
676 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H