11 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
18 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType;
20 template<
typename MatrixType>
21 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
23 typedef typename MatrixType::PlainObject ReturnType;
53 typedef _MatrixType MatrixType;
55 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
57 Options = MatrixType::Options,
58 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
59 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
61 typedef typename MatrixType::Scalar Scalar;
62 typedef typename MatrixType::RealScalar RealScalar;
63 typedef typename MatrixType::Index Index;
64 typedef internal::FullPivHouseholderQRMatrixQReturnType<MatrixType> MatrixQReturnType;
65 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
68 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
69 typedef typename internal::plain_row_type<MatrixType>::type RowVectorType;
70 typedef typename internal::plain_col_type<MatrixType>::type ColVectorType;
80 m_rows_transpositions(),
81 m_cols_transpositions(),
84 m_isInitialized(false),
85 m_usePrescribedThreshold(false) {}
95 m_hCoeffs((std::min)(rows,cols)),
96 m_rows_transpositions(rows),
97 m_cols_transpositions(cols),
98 m_cols_permutation(cols),
99 m_temp((std::min)(rows,cols)),
100 m_isInitialized(false),
101 m_usePrescribedThreshold(false) {}
104 : m_qr(matrix.rows(), matrix.cols()),
105 m_hCoeffs((std::min)(matrix.rows(), matrix.cols())),
106 m_rows_transpositions(matrix.rows()),
107 m_cols_transpositions(matrix.cols()),
108 m_cols_permutation(matrix.cols()),
109 m_temp((std::min)(matrix.rows(), matrix.cols())),
110 m_isInitialized(false),
111 m_usePrescribedThreshold(false)
133 template<
typename Rhs>
134 inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
137 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
138 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*
this, b.derived());
143 MatrixQReturnType
matrixQ(
void)
const;
149 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
155 const PermutationType& colsPermutation()
const
157 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
158 return m_cols_permutation;
161 const IntColVectorType& rowsTranspositions()
const
163 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
164 return m_rows_transpositions;
204 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
207 for(Index i = 0; i < m_nonzero_pivots; ++i)
208 result += (
internal::abs(m_qr.coeff(i,i)) > premultiplied_threshold);
220 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
221 return cols() -
rank();
233 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
234 return rank() == cols();
246 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
247 return rank() == rows();
258 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
267 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
270 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
271 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
272 (*
this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
275 inline Index rows()
const {
return m_qr.rows(); }
276 inline Index cols()
const {
return m_qr.cols(); }
277 const HCoeffsType& hCoeffs()
const {
return m_hCoeffs; }
298 m_usePrescribedThreshold =
true;
313 m_usePrescribedThreshold =
false;
323 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
324 return m_usePrescribedThreshold ? m_prescribedThreshold
339 eigen_assert(m_isInitialized &&
"LU is not initialized.");
340 return m_nonzero_pivots;
350 HCoeffsType m_hCoeffs;
351 IntColVectorType m_rows_transpositions;
352 IntRowVectorType m_cols_transpositions;
353 PermutationType m_cols_permutation;
354 RowVectorType m_temp;
355 bool m_isInitialized, m_usePrescribedThreshold;
356 RealScalar m_prescribedThreshold, m_maxpivot;
357 Index m_nonzero_pivots;
358 RealScalar m_precision;
362 template<
typename MatrixType>
365 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
366 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
370 template<
typename MatrixType>
373 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
374 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
375 return m_qr.diagonal().cwiseAbs().array().log().sum();
378 template<
typename MatrixType>
381 Index rows = matrix.rows();
382 Index cols = matrix.cols();
383 Index size = (std::min)(rows,cols);
386 m_hCoeffs.resize(size);
392 m_rows_transpositions.resize(matrix.rows());
393 m_cols_transpositions.resize(matrix.cols());
394 Index number_of_transpositions = 0;
396 RealScalar biggest(0);
398 m_nonzero_pivots = size;
399 m_maxpivot = RealScalar(0);
401 for (Index k = 0; k < size; ++k)
403 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
404 RealScalar biggest_in_corner;
406 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
408 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
409 row_of_biggest_in_corner += k;
410 col_of_biggest_in_corner += k;
411 if(k==0) biggest = biggest_in_corner;
414 if(internal::isMuchSmallerThan(biggest_in_corner, biggest, m_precision))
416 m_nonzero_pivots = k;
417 for(Index i = k; i < size; i++)
419 m_rows_transpositions.coeffRef(i) = i;
420 m_cols_transpositions.coeffRef(i) = i;
421 m_hCoeffs.coeffRef(i) = Scalar(0);
426 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
427 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
428 if(k != row_of_biggest_in_corner) {
429 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
430 ++number_of_transpositions;
432 if(k != col_of_biggest_in_corner) {
433 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
434 ++number_of_transpositions;
438 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
439 m_qr.coeffRef(k,k) = beta;
444 m_qr.bottomRightCorner(rows-k, cols-k-1)
445 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
448 m_cols_permutation.setIdentity(cols);
449 for(Index k = 0; k < size; ++k)
450 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
452 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
453 m_isInitialized =
true;
460 template<
typename _MatrixType,
typename Rhs>
461 struct solve_retval<FullPivHouseholderQR<_MatrixType>, Rhs>
462 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
464 EIGEN_MAKE_SOLVE_HELPERS(FullPivHouseholderQR<_MatrixType>,Rhs)
466 template<typename Dest>
void evalTo(Dest& dst)
const
468 const Index rows = dec().rows(), cols = dec().cols();
469 eigen_assert(rhs().rows() == rows);
479 typename Rhs::PlainObject c(rhs());
481 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
482 for (Index k = 0; k < dec().rank(); ++k)
484 Index remainingSize = rows-k;
485 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
486 c.bottomRightCorner(remainingSize, rhs().cols())
487 .applyHouseholderOnTheLeft(dec().matrixQR().col(k).tail(remainingSize-1),
488 dec().hCoeffs().coeff(k), &temp.coeffRef(0));
491 if(!dec().isSurjective())
494 RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
495 RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
497 const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
499 if(!internal::isMuchSmallerThan(biggest_in_lower_part_of_c, biggest_in_upper_part_of_c, m_precision))
503 .topLeftCorner(dec().rank(), dec().rank())
504 .template triangularView<Upper>()
505 .solveInPlace(c.topRows(dec().rank()));
507 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
508 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
518 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
519 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
522 typedef typename MatrixType::Index Index;
523 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
524 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
525 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
526 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
528 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
529 const HCoeffsType& hCoeffs,
530 const IntColVectorType& rowsTranspositions)
533 m_rowsTranspositions(rowsTranspositions)
536 template <
typename ResultType>
537 void evalTo(ResultType& result)
const
539 const Index rows = m_qr.rows();
540 WorkVectorType workspace(rows);
541 evalTo(result, workspace);
544 template <
typename ResultType>
545 void evalTo(ResultType& result, WorkVectorType& workspace)
const
550 const Index rows = m_qr.rows();
551 const Index cols = m_qr.cols();
552 const Index size = (std::min)(rows, cols);
553 workspace.resize(rows);
554 result.setIdentity(rows, rows);
555 for (Index k = size-1; k >= 0; k--)
557 result.block(k, k, rows-k, rows-k)
558 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
559 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
563 Index rows()
const {
return m_qr.rows(); }
564 Index cols()
const {
return m_qr.rows(); }
567 typename MatrixType::Nested m_qr;
568 typename HCoeffsType::Nested m_hCoeffs;
569 typename IntColVectorType::Nested m_rowsTranspositions;
574 template<
typename MatrixType>
577 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
578 return MatrixQReturnType(m_qr, m_hCoeffs, m_rows_transpositions);
585 template<
typename Derived>
594 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H