26 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
27 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H
33 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType;
35 template<
typename MatrixType>
36 struct traits<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
38 typedef typename MatrixType::PlainObject ReturnType;
76 typedef typename MatrixType::Scalar
Scalar;
78 typedef typename MatrixType::Index
Index;
80 typedef typename internal::plain_diag_type<MatrixType>::type
HCoeffsType;
84 typedef typename internal::plain_row_type<MatrixType>::type
RowVectorType;
85 typedef typename internal::plain_col_type<MatrixType>::type
ColVectorType;
114 m_temp((std::min)(rows,cols)),
148 template<
typename Rhs>
149 inline const internal::solve_retval<FullPivHouseholderQR, Rhs>
153 return internal::solve_retval<FullPivHouseholderQR, Rhs>(*
this, b.derived());
282 internal::solve_retval<FullPivHouseholderQR, typename MatrixType::IdentityReturnType>
286 return internal::solve_retval<FullPivHouseholderQR,typename MatrixType::IdentityReturnType>
287 (*
this, MatrixType::Identity(
m_qr.rows(),
m_qr.cols()));
377 template<
typename MatrixType>
380 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
381 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
385 template<
typename MatrixType>
388 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
389 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
390 return m_qr.diagonal().cwiseAbs().array().log().sum();
393 template<
typename MatrixType>
397 Index cols = matrix.cols();
398 Index size = (std::min)(rows,cols);
401 m_hCoeffs.resize(size);
407 m_rows_transpositions.resize(matrix.rows());
408 m_cols_transpositions.resize(matrix.cols());
409 Index number_of_transpositions = 0;
413 m_nonzero_pivots = size;
416 for (
Index k = 0; k < size; ++k)
418 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
421 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
423 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
424 row_of_biggest_in_corner += k;
425 col_of_biggest_in_corner += k;
426 if(k==0) biggest = biggest_in_corner;
431 m_nonzero_pivots = k;
432 for(
Index i = k; i < size; i++)
434 m_rows_transpositions.coeffRef(i) = i;
435 m_cols_transpositions.coeffRef(i) = i;
436 m_hCoeffs.coeffRef(i) =
Scalar(0);
441 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
442 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
443 if(k != row_of_biggest_in_corner) {
444 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
445 ++number_of_transpositions;
447 if(k != col_of_biggest_in_corner) {
448 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
449 ++number_of_transpositions;
453 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
454 m_qr.coeffRef(k,k) = beta;
459 m_qr.bottomRightCorner(rows-k, cols-k-1)
460 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
463 m_cols_permutation.setIdentity(cols);
464 for(
Index k = 0; k < size; ++k)
465 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
467 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
468 m_isInitialized =
true;
475 template<
typename _MatrixType,
typename Rhs>
477 : solve_retval_base<FullPivHouseholderQR<_MatrixType>, Rhs>
481 template<typename Dest>
void evalTo(Dest& dst)
const
483 const Index rows = dec().rows(), cols = dec().cols();
494 typename Rhs::PlainObject c(rhs());
496 Matrix<Scalar,1,Rhs::ColsAtCompileTime> temp(rhs().cols());
497 for (Index k = 0; k < dec().rank(); ++k)
499 Index remainingSize = rows-k;
500 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
501 c.bottomRightCorner(remainingSize, rhs().cols())
502 .applyHouseholderOnTheLeft(dec().matrixQR().
col(k).tail(remainingSize-1),
503 dec().hCoeffs().coeff(k), &temp.coeffRef(0));
506 if(!dec().isSurjective())
509 RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
510 RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
512 const RealScalar m_precision = NumTraits<Scalar>::epsilon() * (std::min)(rows,cols);
518 .topLeftCorner(dec().rank(), dec().rank())
519 .template triangularView<Upper>()
520 .solveInPlace(c.topRows(dec().rank()));
522 for(Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
523 for(Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
533 template<
typename MatrixType>
struct FullPivHouseholderQRMatrixQReturnType
534 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
537 typedef typename MatrixType::Index Index;
538 typedef typename internal::plain_col_type<MatrixType, Index>::type IntColVectorType;
539 typedef typename internal::plain_diag_type<MatrixType>::type HCoeffsType;
540 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
541 MatrixType::MaxRowsAtCompileTime> WorkVectorType;
543 FullPivHouseholderQRMatrixQReturnType(
const MatrixType& qr,
544 const HCoeffsType& hCoeffs,
545 const IntColVectorType& rowsTranspositions)
548 m_rowsTranspositions(rowsTranspositions)
551 template <
typename ResultType>
552 void evalTo(ResultType& result)
const
554 const Index rows = m_qr.rows();
555 WorkVectorType workspace(rows);
556 evalTo(result, workspace);
559 template <
typename ResultType>
560 void evalTo(ResultType& result, WorkVectorType& workspace)
const
565 const Index rows = m_qr.rows();
566 const Index cols = m_qr.cols();
567 const Index size = (std::min)(rows, cols);
568 workspace.resize(rows);
569 result.setIdentity(rows, rows);
570 for (Index k = size-1; k >= 0; k--)
572 result.block(k, k, rows-k, rows-k)
573 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
internal::conj(m_hCoeffs.coeff(k)), &workspace.coeffRef(k));
574 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
578 Index rows()
const {
return m_qr.rows(); }
579 Index cols()
const {
return m_qr.rows(); }
582 typename MatrixType::Nested m_qr;
583 typename HCoeffsType::Nested m_hCoeffs;
584 typename IntColVectorType::Nested m_rowsTranspositions;
589 template<
typename MatrixType>
592 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
600 template<
typename Derived>
609 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H