11 #ifndef EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 12 #define EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H 20 template<
typename MatrixType>
55 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
56 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
58 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
59 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
61 typedef typename MatrixType::Scalar
Scalar;
63 typedef typename MatrixType::Index
Index;
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) {}
116 : m_qr(matrix.rows(), matrix.cols()),
117 m_hCoeffs((
std::min)(matrix.rows(), matrix.cols())),
118 m_rows_transpositions(matrix.rows()),
119 m_cols_transpositions(matrix.cols()),
120 m_cols_permutation(matrix.cols()),
121 m_temp((
std::min)(matrix.rows(), matrix.cols())),
122 m_isInitialized(false),
123 m_usePrescribedThreshold(false)
145 template<
typename Rhs>
149 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
155 MatrixQReturnType matrixQ(
void)
const;
161 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
170 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
171 return m_cols_permutation;
177 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
178 return m_rows_transpositions;
194 typename MatrixType::RealScalar absDeterminant()
const;
208 typename MatrixType::RealScalar logAbsDeterminant()
const;
219 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
220 RealScalar premultiplied_threshold =
abs(m_maxpivot) * threshold();
222 for(Index i = 0; i < m_nonzero_pivots; ++i)
223 result += (
abs(m_qr.coeff(i,i)) > premultiplied_threshold);
235 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
236 return cols() - rank();
248 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
249 return rank() == cols();
261 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
262 return rank() == rows();
273 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
274 return isInjective() && isSurjective();
285 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
287 (*
this, MatrixType::Identity(m_qr.rows(), m_qr.cols()));
290 inline Index
rows()
const {
return m_qr.rows(); }
291 inline Index
cols()
const {
return m_qr.cols(); }
297 const HCoeffsType&
hCoeffs()
const {
return m_hCoeffs; }
318 m_usePrescribedThreshold =
true;
319 m_prescribedThreshold = threshold;
333 m_usePrescribedThreshold =
false;
343 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
344 return m_usePrescribedThreshold ? m_prescribedThreshold
359 eigen_assert(m_isInitialized &&
"LU is not initialized.");
360 return m_nonzero_pivots;
382 template<
typename MatrixType>
386 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
387 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
388 return abs(m_qr.diagonal().prod());
391 template<
typename MatrixType>
394 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
395 eigen_assert(m_qr.rows() == m_qr.cols() &&
"You can't take the determinant of a non-square matrix!");
396 return m_qr.diagonal().cwiseAbs().array().log().sum();
405 template<
typename MatrixType>
409 Index rows = matrix.rows();
410 Index cols = matrix.cols();
411 Index size = (std::min)(rows,cols);
414 m_hCoeffs.resize(size);
420 m_rows_transpositions.resize(matrix.rows());
421 m_cols_transpositions.resize(matrix.cols());
422 Index number_of_transpositions = 0;
426 m_nonzero_pivots = size;
429 for (
Index k = 0; k < size; ++k)
431 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
434 biggest_in_corner = m_qr.bottomRightCorner(rows-k, cols-k)
436 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
437 row_of_biggest_in_corner += k;
438 col_of_biggest_in_corner += k;
439 if(k==0) biggest = biggest_in_corner;
444 m_nonzero_pivots = k;
445 for(
Index i = k; i < size; i++)
447 m_rows_transpositions.coeffRef(i) = i;
448 m_cols_transpositions.coeffRef(i) = i;
449 m_hCoeffs.coeffRef(i) =
Scalar(0);
454 m_rows_transpositions.coeffRef(k) = row_of_biggest_in_corner;
455 m_cols_transpositions.coeffRef(k) = col_of_biggest_in_corner;
456 if(k != row_of_biggest_in_corner) {
457 m_qr.row(k).tail(cols-k).swap(m_qr.row(row_of_biggest_in_corner).tail(cols-k));
458 ++number_of_transpositions;
460 if(k != col_of_biggest_in_corner) {
461 m_qr.col(k).swap(m_qr.col(col_of_biggest_in_corner));
462 ++number_of_transpositions;
466 m_qr.col(k).tail(rows-k).makeHouseholderInPlace(m_hCoeffs.coeffRef(k), beta);
467 m_qr.coeffRef(k,k) = beta;
470 if(
abs(beta) > m_maxpivot) m_maxpivot =
abs(beta);
472 m_qr.bottomRightCorner(rows-k, cols-k-1)
473 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1), m_hCoeffs.coeffRef(k), &m_temp.coeffRef(k+1));
476 m_cols_permutation.setIdentity(cols);
477 for(
Index k = 0; k < size; ++k)
478 m_cols_permutation.applyTranspositionOnTheRight(k, m_cols_transpositions.coeff(k));
480 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
481 m_isInitialized =
true;
488 template<
typename _MatrixType,
typename Rhs>
494 template<typename Dest>
void evalTo(Dest& dst)
const 496 const Index rows = dec().rows(), cols = dec().cols();
507 typename Rhs::PlainObject c(
rhs());
510 for (
Index k = 0; k < dec().rank(); ++k)
512 Index remainingSize = rows-k;
513 c.row(k).swap(c.row(dec().rowsTranspositions().coeff(k)));
514 c.bottomRightCorner(remainingSize,
rhs().cols())
515 .applyHouseholderOnTheLeft(dec().matrixQR().
col(k).
tail(remainingSize-1),
516 dec().hCoeffs().coeff(k), &temp.
coeffRef(0));
519 if(!dec().isSurjective())
522 RealScalar biggest_in_upper_part_of_c = c.topRows( dec().rank() ).cwiseAbs().maxCoeff();
523 RealScalar biggest_in_lower_part_of_c = c.bottomRows(rows-dec().rank()).cwiseAbs().maxCoeff();
531 .topLeftCorner(dec().rank(), dec().rank())
532 .template triangularView<Upper>()
533 .solveInPlace(c.topRows(dec().rank()));
535 for(
Index i = 0; i < dec().rank(); ++i) dst.row(dec().colsPermutation().indices().coeff(i)) = c.row(i);
536 for(
Index i = dec().rank(); i < cols; ++i) dst.row(dec().colsPermutation().indices().coeff(i)).setZero();
547 :
public ReturnByValue<FullPivHouseholderQRMatrixQReturnType<MatrixType> >
550 typedef typename MatrixType::Index
Index;
553 typedef Matrix<
typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime,
RowMajor, 1,
557 const HCoeffsType& hCoeffs,
558 const IntColVectorType& rowsTranspositions)
561 m_rowsTranspositions(rowsTranspositions)
564 template <
typename ResultType>
567 const Index rows = m_qr.rows();
568 WorkVectorType workspace(rows);
569 evalTo(result, workspace);
572 template <
typename ResultType>
573 void evalTo(ResultType& result, WorkVectorType& workspace)
const 579 const Index rows = m_qr.rows();
580 const Index cols = m_qr.cols();
581 const Index size = (std::min)(rows, cols);
583 result.setIdentity(rows, rows);
584 for (Index k = size-1; k >= 0; k--)
586 result.block(k, k, rows-k, rows-k)
587 .applyHouseholderOnTheLeft(m_qr.col(k).tail(rows-k-1),
conj(m_hCoeffs.coeff(k)), &workspace.
coeffRef(k));
588 result.row(k).swap(result.row(m_rowsTranspositions.coeff(k)));
592 Index
rows()
const {
return m_qr.rows(); }
593 Index
cols()
const {
return m_qr.rows(); }
596 typename MatrixType::Nested
m_qr;
603 template<
typename MatrixType>
606 eigen_assert(m_isInitialized &&
"FullPivHouseholderQR is not initialized.");
614 template<
typename Derived>
623 #endif // EIGEN_FULLPIVOTINGHOUSEHOLDERQR_H Index dimensionOfKernel() const
bool isInvertible() const
FullPivHouseholderQR & setThreshold(Default_t)
const AutoDiffScalar< DerType > & conj(const AutoDiffScalar< DerType > &x)
Householder rank-revealing QR decomposition of a matrix with full pivoting.
internal::plain_col_type< MatrixType, Index >::type IntColVectorType
const internal::solve_retval< FullPivHouseholderQR, typename MatrixType::IdentityReturnType > inverse() const
void evalTo(ResultType &result) const
const IntColVectorType & rowsTranspositions() const
IntColVectorType m_rows_transpositions
MatrixQReturnType matrixQ(void) const
bool m_usePrescribedThreshold
RealScalar maxPivot() const
iterative scaling algorithm to equilibrate rows and column norms in matrices
const internal::solve_retval< FullPivHouseholderQR, Rhs > solve(const MatrixBase< Rhs > &b) const
bool isSurjective() const
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
MatrixType::RealScalar RealScalar
bool isMuchSmallerThan(const Scalar &x, const OtherScalar &y, typename NumTraits< Scalar >::Real precision=NumTraits< Scalar >::dummy_precision())
void evalTo(ResultType &result, WorkVectorType &workspace) const
FullPivHouseholderQRMatrixQReturnType(const MatrixType &qr, const HCoeffsType &hCoeffs, const IntColVectorType &rowsTranspositions)
IntRowVectorType m_cols_transpositions
FullPivHouseholderQR(Index rows, Index cols)
Default Constructor with memory preallocation.
internal::plain_diag_type< MatrixType >::type HCoeffsType
EIGEN_STRONG_INLINE const CwiseUnaryOp< internal::scalar_abs_op< Scalar >, const Derived > abs() const
Matrix< typename MatrixType::Scalar, 1, MatrixType::RowsAtCompileTime, RowMajor, 1, MatrixType::MaxRowsAtCompileTime > WorkVectorType
const HCoeffsType & hCoeffs() const
HCoeffsType::Nested m_hCoeffs
FullPivHouseholderQR()
Default Constructor.
internal::plain_row_type< MatrixType >::type RowVectorType
const PermutationType & colsPermutation() const
internal::plain_col_type< MatrixType, Index >::type IntColVectorType
FullPivHouseholderQR(const MatrixType &matrix)
Constructs a QR factorization from a given matrix.
PermutationType m_cols_permutation
EIGEN_STRONG_INLINE Scalar & coeffRef(Index rowId, Index colId)
EIGEN_STRONG_INLINE void resize(Index nbRows, Index nbCols)
SegmentReturnType tail(Index vecSize)
Provides a generic way to set and pass user-specified options.
internal::FullPivHouseholderQRMatrixQReturnType< MatrixType > MatrixQReturnType
Expression type for return value of FullPivHouseholderQR::matrixQ()
internal::plain_col_type< MatrixType >::type ColVectorType
FullPivHouseholderQR & setThreshold(const RealScalar &threshold)
MatrixType::Scalar Scalar
PermutationMatrix< ColsAtCompileTime, MaxColsAtCompileTime > PermutationType
void rhs(const real_t *x, real_t *f)
MatrixType::RealScalar logAbsDeterminant() const
RealScalar threshold() const
MatrixType::RealScalar absDeterminant() const
const MatrixType & matrixQR() const
Matrix< Index, 1, ColsAtCompileTime, RowMajor, 1, MaxColsAtCompileTime > IntRowVectorType
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
RealScalar m_prescribedThreshold
MatrixType::PlainObject ReturnType
FullPivHouseholderQR & compute(const MatrixType &matrix)
Base class for all dense matrices, vectors, and expressions.
internal::plain_diag_type< MatrixType >::type HCoeffsType
const FullPivHouseholderQR< PlainObject > fullPivHouseholderQr() const
Index nonzeroPivots() const
IntColVectorType::Nested m_rowsTranspositions