16 template<
typename _MatrixType>
struct traits<FullPivLU<_MatrixType> >
19 typedef MatrixXpr XprKind;
20 typedef SolverStorage StorageKind;
57 template<
typename _MatrixType>
class FullPivLU
58 :
public SolverBase<FullPivLU<_MatrixType> >
61 typedef _MatrixType MatrixType;
62 typedef SolverBase<FullPivLU> Base;
64 EIGEN_GENERIC_PUBLIC_INTERFACE(FullPivLU)
67 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
68 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
70 typedef typename internal::plain_row_type<MatrixType, StorageIndex>::type IntRowVectorType;
71 typedef typename internal::plain_col_type<MatrixType, StorageIndex>::type IntColVectorType;
72 typedef PermutationMatrix<ColsAtCompileTime, MaxColsAtCompileTime> PermutationQType;
73 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationPType;
74 typedef typename MatrixType::PlainObject PlainObject;
90 FullPivLU(Index rows, Index cols);
97 template<
typename InputType>
98 explicit FullPivLU(
const EigenBase<InputType>& matrix);
107 template<
typename InputType>
108 FullPivLU& compute(
const EigenBase<InputType>& matrix);
118 eigen_assert(m_isInitialized &&
"LU is not initialized.");
131 eigen_assert(m_isInitialized &&
"LU is not initialized.");
132 return m_nonzero_pivots;
146 eigen_assert(m_isInitialized &&
"LU is not initialized.");
156 eigen_assert(m_isInitialized &&
"LU is not initialized.");
174 inline const internal::kernel_retval<FullPivLU>
kernel()
const 176 eigen_assert(m_isInitialized &&
"LU is not initialized.");
177 return internal::kernel_retval<FullPivLU>(*this);
199 inline const internal::image_retval<FullPivLU>
200 image(
const MatrixType& originalMatrix)
const 202 eigen_assert(m_isInitialized &&
"LU is not initialized.");
203 return internal::image_retval<FullPivLU>(*
this, originalMatrix);
226 template<
typename Rhs>
230 eigen_assert(m_isInitialized &&
"LU is not initialized.");
249 typename internal::traits<MatrixType>::Scalar determinant()
const;
270 m_usePrescribedThreshold =
true;
271 m_prescribedThreshold = threshold;
285 m_usePrescribedThreshold =
false;
295 eigen_assert(m_isInitialized || m_usePrescribedThreshold);
296 return m_usePrescribedThreshold ? m_prescribedThreshold
311 eigen_assert(m_isInitialized &&
"LU is not initialized.");
312 RealScalar premultiplied_threshold = abs(m_maxpivot) * threshold();
314 for(
Index i = 0; i < m_nonzero_pivots; ++i)
315 result += (abs(m_lu.coeff(i,i)) > premultiplied_threshold);
327 eigen_assert(m_isInitialized &&
"LU is not initialized.");
328 return cols() - rank();
340 eigen_assert(m_isInitialized &&
"LU is not initialized.");
341 return rank() == cols();
353 eigen_assert(m_isInitialized &&
"LU is not initialized.");
354 return rank() == rows();
365 eigen_assert(m_isInitialized &&
"LU is not initialized.");
366 return isInjective() && (m_lu.rows() == m_lu.cols());
378 eigen_assert(m_isInitialized &&
"LU is not initialized.");
379 eigen_assert(m_lu.rows() == m_lu.cols() &&
"You can't take the inverse of a non-square matrix!");
383 MatrixType reconstructedMatrix()
const;
385 inline Index rows()
const {
return m_lu.rows(); }
386 inline Index cols()
const {
return m_lu.cols(); }
388 #ifndef EIGEN_PARSED_BY_DOXYGEN 389 template<
typename RhsType,
typename DstType>
391 void _solve_impl(
const RhsType &rhs, DstType &dst)
const;
393 template<
bool Conjugate,
typename RhsType,
typename DstType>
395 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const;
400 static void check_template_parameters()
402 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
405 void computeInPlace();
408 PermutationPType m_p;
409 PermutationQType m_q;
410 IntColVectorType m_rowsTranspositions;
411 IntRowVectorType m_colsTranspositions;
412 Index m_det_pq, m_nonzero_pivots;
413 RealScalar m_maxpivot, m_prescribedThreshold;
414 bool m_isInitialized, m_usePrescribedThreshold;
417 template<
typename MatrixType>
419 : m_isInitialized(false), m_usePrescribedThreshold(false)
423 template<
typename MatrixType>
428 m_rowsTranspositions(rows),
429 m_colsTranspositions(cols),
430 m_isInitialized(false),
431 m_usePrescribedThreshold(false)
435 template<
typename MatrixType>
436 template<
typename InputType>
438 : m_lu(matrix.rows(), matrix.cols()),
441 m_rowsTranspositions(matrix.rows()),
442 m_colsTranspositions(matrix.cols()),
443 m_isInitialized(false),
444 m_usePrescribedThreshold(false)
449 template<
typename MatrixType>
450 template<
typename InputType>
453 check_template_parameters();
458 m_isInitialized =
true;
466 template<
typename MatrixType>
470 const Index rows = m_lu.rows();
471 const Index cols = m_lu.cols();
475 m_rowsTranspositions.resize(m_lu.rows());
476 m_colsTranspositions.resize(m_lu.cols());
477 Index number_of_transpositions = 0;
479 m_nonzero_pivots =
size;
480 m_maxpivot = RealScalar(0);
487 Index row_of_biggest_in_corner, col_of_biggest_in_corner;
488 typedef internal::scalar_score_coeff_op<Scalar> Scoring;
489 typedef typename Scoring::result_type Score;
490 Score biggest_in_corner;
491 biggest_in_corner = m_lu.bottomRightCorner(rows-k, cols-k)
492 .unaryExpr(Scoring())
493 .maxCoeff(&row_of_biggest_in_corner, &col_of_biggest_in_corner);
494 row_of_biggest_in_corner += k;
495 col_of_biggest_in_corner += k;
497 if(biggest_in_corner==Score(0))
501 m_nonzero_pivots = k;
504 m_rowsTranspositions.coeffRef(i) = i;
505 m_colsTranspositions.coeffRef(i) = i;
510 RealScalar abs_pivot = internal::abs_knowing_score<Scalar>()(m_lu(row_of_biggest_in_corner, col_of_biggest_in_corner), biggest_in_corner);
511 if(abs_pivot > m_maxpivot) m_maxpivot = abs_pivot;
516 m_rowsTranspositions.coeffRef(k) = row_of_biggest_in_corner;
517 m_colsTranspositions.coeffRef(k) = col_of_biggest_in_corner;
518 if(k != row_of_biggest_in_corner) {
519 m_lu.row(k).swap(m_lu.row(row_of_biggest_in_corner));
520 ++number_of_transpositions;
522 if(k != col_of_biggest_in_corner) {
523 m_lu.col(k).swap(m_lu.col(col_of_biggest_in_corner));
524 ++number_of_transpositions;
531 m_lu.col(k).tail(rows-k-1) /= m_lu.coeff(k,k);
533 m_lu.block(k+1,k+1,rows-k-1,cols-k-1).noalias() -= m_lu.col(k).tail(rows-k-1) * m_lu.row(k).tail(cols-k-1);
540 for(
Index k = size-1; k >= 0; --k)
547 m_det_pq = (number_of_transpositions%2) ? -1 : 1;
550 template<
typename MatrixType>
553 eigen_assert(m_isInitialized &&
"LU is not initialized.");
554 eigen_assert(m_lu.rows() == m_lu.cols() &&
"You can't take the determinant of a non-square matrix!");
555 return Scalar(m_det_pq) * Scalar(m_lu.diagonal().prod());
561 template<
typename MatrixType>
564 eigen_assert(m_isInitialized &&
"LU is not initialized.");
565 const Index smalldim = (std::min)(m_lu.rows(), m_lu.cols());
567 MatrixType res(m_lu.rows(),m_lu.cols());
569 res = m_lu.leftCols(smalldim)
570 .template triangularView<UnitLower>().toDenseMatrix()
571 * m_lu.topRows(smalldim)
572 .template triangularView<Upper>().toDenseMatrix();
586 template<
typename _MatrixType>
587 struct kernel_retval<FullPivLU<_MatrixType> >
588 : kernel_retval_base<FullPivLU<_MatrixType> >
592 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
593 MatrixType::MaxColsAtCompileTime,
594 MatrixType::MaxRowsAtCompileTime)
597 template<
typename Dest>
void evalTo(Dest& dst)
const 600 const Index cols = dec().matrixLU().cols(), dimker = cols -
rank();
627 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
629 for(
Index i = 0; i < dec().nonzeroPivots(); ++i)
630 if(abs(dec().
matrixLU().coeff(i,i)) > premultiplied_threshold)
631 pivots.coeffRef(p++) = i;
632 eigen_internal_assert(p ==
rank());
638 Matrix<
typename MatrixType::Scalar, Dynamic, Dynamic, MatrixType::Options,
639 MaxSmallDimAtCompileTime, MatrixType::MaxColsAtCompileTime>
643 if(i) m.row(i).head(i).setZero();
644 m.row(i).tail(cols-i) = dec().matrixLU().row(pivots.coeff(i)).tail(cols-i);
647 m.block(0, 0,
rank(),
rank()).template triangularView<StrictlyLower>().setZero();
649 m.col(i).swap(m.col(pivots.coeff(i)));
655 .template triangularView<Upper>().solveInPlace(
656 m.topRightCorner(rank(), dimker)
661 m.col(i).swap(m.col(pivots.coeff(i)));
664 for(
Index i = 0; i <
rank(); ++i) dst.row(dec().
permutationQ().indices().coeff(i)) = -m.row(i).tail(dimker);
666 for(
Index k = 0; k < dimker; ++k) dst.coeffRef(dec().
permutationQ().indices().coeff(
rank()+k), k) = Scalar(1);
672 template<
typename _MatrixType>
673 struct image_retval<FullPivLU<_MatrixType> >
674 : image_retval_base<FullPivLU<_MatrixType> >
678 enum { MaxSmallDimAtCompileTime = EIGEN_SIZE_MIN_PREFER_FIXED(
679 MatrixType::MaxColsAtCompileTime,
680 MatrixType::MaxRowsAtCompileTime)
683 template<
typename Dest>
void evalTo(Dest& dst)
const 696 RealScalar premultiplied_threshold = dec().maxPivot() * dec().threshold();
698 for(
Index i = 0; i < dec().nonzeroPivots(); ++i)
699 if(abs(dec().
matrixLU().coeff(i,i)) > premultiplied_threshold)
700 pivots.coeffRef(p++) = i;
701 eigen_internal_assert(p ==
rank());
704 dst.col(i) = originalMatrix().col(dec().
permutationQ().indices().coeff(pivots.coeff(i)));
712 #ifndef EIGEN_PARSED_BY_DOXYGEN 713 template<
typename _MatrixType>
714 template<
typename RhsType,
typename DstType>
725 const Index rows = this->rows(),
727 nonzero_pivots = this->
rank();
728 eigen_assert(rhs.rows() == rows);
729 const Index smalldim = (std::min)(rows, cols);
731 if(nonzero_pivots == 0)
737 typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
743 m_lu.topLeftCorner(smalldim,smalldim)
744 .template triangularView<UnitLower>()
745 .solveInPlace(c.topRows(smalldim));
747 c.bottomRows(rows-cols) -= m_lu.bottomRows(rows-cols) * c.topRows(cols);
750 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
751 .template triangularView<Upper>()
752 .solveInPlace(c.topRows(nonzero_pivots));
755 for(
Index i = 0; i < nonzero_pivots; ++i)
757 for(
Index i = nonzero_pivots; i < m_lu.cols(); ++i)
761 template<
typename _MatrixType>
762 template<
bool Conjugate,
typename RhsType,
typename DstType>
776 const Index rows = this->rows(), cols = this->cols(),
777 nonzero_pivots = this->
rank();
778 eigen_assert(rhs.rows() == cols);
779 const Index smalldim = (std::min)(rows, cols);
781 if(nonzero_pivots == 0)
787 typename RhsType::PlainObject c(rhs.rows(), rhs.cols());
794 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
795 .template triangularView<Upper>()
797 .solveInPlace(c.topRows(nonzero_pivots));
799 m_lu.topLeftCorner(smalldim, smalldim)
800 .template triangularView<UnitLower>()
802 .solveInPlace(c.topRows(smalldim));
805 m_lu.topLeftCorner(nonzero_pivots, nonzero_pivots)
806 .template triangularView<Upper>()
808 .solveInPlace(c.topRows(nonzero_pivots));
810 m_lu.topLeftCorner(smalldim, smalldim)
811 .template triangularView<UnitLower>()
813 .solveInPlace(c.topRows(smalldim));
818 for(
Index i = 0; i < smalldim; ++i)
819 dst.row(invp.
indices().coeff(i)) = c.row(i);
820 for(
Index i = smalldim; i < rows; ++i)
821 dst.row(invp.
indices().coeff(i)).setZero();
830 template<
typename DstXprType,
typename MatrixType,
typename Scalar>
831 struct Assignment<DstXprType, Inverse<FullPivLU<MatrixType> >, internal::assign_op<Scalar>, Dense2Dense, Scalar>
835 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
837 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
851 template<
typename Derived>
bool isInjective() const
Definition: FullPivLU.h:338
const internal::image_retval< FullPivLU > image(const MatrixType &originalMatrix) const
Definition: FullPivLU.h:200
FullPivLU & setThreshold(const RealScalar &threshold)
Definition: FullPivLU.h:268
internal::traits< MatrixType >::Scalar determinant() const
Definition: FullPivLU.h:551
FullPivLU & compute(const EigenBase< InputType > &matrix)
Derived & applyTranspositionOnTheRight(Index i, Index j)
Definition: PermutationMatrix.h:186
AdjointReturnType adjoint() const
Definition: SolverBase.h:109
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:107
Derived & derived()
Definition: EigenBase.h:44
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:37
Index rows() const
Definition: EigenBase.h:58
RealScalar threshold() const
Definition: FullPivLU.h:293
Index nonzeroPivots() const
Definition: FullPivLU.h:129
Definition: EigenBase.h:28
Expression of the inverse of another expression.
Definition: Inverse.h:43
bool isSurjective() const
Definition: FullPivLU.h:351
FullPivLU & setThreshold(Default_t)
Definition: FullPivLU.h:283
const PermutationQType & permutationQ() const
Definition: FullPivLU.h:154
void setIdentity()
Definition: PermutationMatrix.h:143
MatrixType reconstructedMatrix() const
Definition: FullPivLU.h:562
ConstTransposeReturnType transpose() const
Definition: SolverBase.h:90
const PermutationPType & permutationP() const
Definition: FullPivLU.h:144
InverseReturnType inverse() const
Definition: PermutationMatrix.h:197
Definition: Eigen_Colamd.h:54
Index cols() const
Definition: EigenBase.h:61
LU decomposition of a matrix with complete pivoting, and related features.
Definition: ForwardDeclarations.h:247
bool isInvertible() const
Definition: FullPivLU.h:363
const FullPivLU< PlainObject > fullPivLu() const
Definition: FullPivLU.h:853
const IndicesType & indices() const
Definition: PermutationMatrix.h:390
FullPivLU()
Default Constructor.
Definition: FullPivLU.h:418
Index rank() const
Definition: FullPivLU.h:308
const MatrixType & matrixLU() const
Definition: FullPivLU.h:116
Pseudo expression representing a solving operation.
Definition: Solve.h:62
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:178
const Inverse< FullPivLU > inverse() const
Definition: FullPivLU.h:376
Index dimensionOfKernel() const
Definition: FullPivLU.h:325
RealScalar maxPivot() const
Definition: FullPivLU.h:138
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const Solve< FullPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: FullPivLU.h:228
Index size() const
Definition: EigenBase.h:65
const internal::kernel_retval< FullPivLU > kernel() const
Definition: FullPivLU.h:174