11 #ifndef EIGEN_PARTIALLU_H 12 #define EIGEN_PARTIALLU_H 17 template<
typename _MatrixType>
struct traits<PartialPivLU<_MatrixType> >
20 typedef MatrixXpr XprKind;
21 typedef SolverStorage StorageKind;
22 typedef traits<_MatrixType> BaseTraits;
29 template<
typename T,
typename Derived>
35 template<
typename T,
typename Derived>
36 struct enable_if_ref<Ref<T>,Derived> {
75 template<
typename _MatrixType>
class PartialPivLU
76 :
public SolverBase<PartialPivLU<_MatrixType> >
80 typedef _MatrixType MatrixType;
81 typedef SolverBase<PartialPivLU> Base;
82 EIGEN_GENERIC_PUBLIC_INTERFACE(PartialPivLU)
85 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
86 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
88 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
89 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
90 typedef typename MatrixType::PlainObject PlainObject;
106 explicit PartialPivLU(
Index size);
115 template<
typename InputType>
116 explicit PartialPivLU(
const EigenBase<InputType>& matrix);
125 template<
typename InputType>
126 explicit PartialPivLU(EigenBase<InputType>& matrix);
128 template<
typename InputType>
129 PartialPivLU& compute(
const EigenBase<InputType>& matrix) {
130 m_lu = matrix.derived();
143 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
151 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
173 template<
typename Rhs>
177 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
186 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
187 return internal::rcond_estimate_helper(m_l1_norm, *
this);
199 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
216 Scalar determinant()
const;
218 MatrixType reconstructedMatrix()
const;
220 inline Index rows()
const {
return m_lu.rows(); }
221 inline Index cols()
const {
return m_lu.cols(); }
223 #ifndef EIGEN_PARSED_BY_DOXYGEN 224 template<
typename RhsType,
typename DstType>
226 void _solve_impl(
const RhsType &rhs, DstType &dst)
const {
234 eigen_assert(rhs.rows() == m_lu.rows());
237 dst = permutationP() * rhs;
240 m_lu.template triangularView<UnitLower>().solveInPlace(dst);
243 m_lu.template triangularView<Upper>().solveInPlace(dst);
246 template<
bool Conjugate,
typename RhsType,
typename DstType>
248 void _solve_impl_transposed(
const RhsType &rhs, DstType &dst)
const {
256 eigen_assert(rhs.rows() == m_lu.cols());
260 dst = m_lu.template triangularView<Upper>().adjoint().solve(rhs);
262 m_lu.template triangularView<UnitLower>().adjoint().solveInPlace(dst);
265 dst = m_lu.template triangularView<Upper>().transpose().solve(rhs);
267 m_lu.template triangularView<UnitLower>().transpose().solveInPlace(dst);
270 dst = permutationP().transpose() * dst;
276 static void check_template_parameters()
278 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
285 TranspositionType m_rowsTranspositions;
286 RealScalar m_l1_norm;
288 bool m_isInitialized;
291 template<
typename MatrixType>
295 m_rowsTranspositions(),
298 m_isInitialized(false)
302 template<
typename MatrixType>
306 m_rowsTranspositions(size),
309 m_isInitialized(false)
313 template<
typename MatrixType>
314 template<
typename InputType>
316 : m_lu(matrix.rows(),matrix.cols()),
318 m_rowsTranspositions(matrix.rows()),
321 m_isInitialized(false)
326 template<
typename MatrixType>
327 template<
typename InputType>
331 m_rowsTranspositions(matrix.rows()),
334 m_isInitialized(false)
342 template<
typename Scalar,
int StorageOrder,
typename PivIndex>
343 struct partial_lu_impl
353 typedef typename MatrixType::RealScalar RealScalar;
365 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
367 typedef scalar_score_coeff_op<Scalar> Scoring;
368 typedef typename Scoring::result_type Score;
369 const Index rows = lu.rows();
370 const Index cols = lu.cols();
371 const Index size = (std::min)(rows,cols);
372 nb_transpositions = 0;
373 Index first_zero_pivot = -1;
376 Index rrows = rows-k-1;
377 Index rcols = cols-k-1;
379 Index row_of_biggest_in_col;
380 Score biggest_in_corner
381 = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
382 row_of_biggest_in_col += k;
384 row_transpositions[k] = PivIndex(row_of_biggest_in_col);
386 if(biggest_in_corner != Score(0))
388 if(k != row_of_biggest_in_col)
390 lu.row(k).swap(lu.row(row_of_biggest_in_col));
396 lu.col(k).tail(rrows) /= lu.coeff(k,k);
398 else if(first_zero_pivot==-1)
402 first_zero_pivot = k;
406 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
408 return first_zero_pivot;
426 static Index blocked_lu(
Index rows,
Index cols, Scalar* lu_data,
Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions,
Index maxBlockSize=256)
428 MapLU lu1(lu_data,StorageOrder==
RowMajor?rows:luStride,StorageOrder==
RowMajor?luStride:cols);
429 MatrixType lu(lu1,0,0,rows,cols);
431 const Index size = (std::min)(rows,cols);
436 return unblocked_lu(lu, row_transpositions, nb_transpositions);
444 blockSize = (blockSize/16)*16;
445 blockSize = (std::min)((std::max)(blockSize,
Index(8)), maxBlockSize);
448 nb_transpositions = 0;
449 Index first_zero_pivot = -1;
452 Index bs = (std::min)(size-k,blockSize);
453 Index trows = rows - k - bs;
454 Index tsize = size - k - bs;
460 BlockType A_0(lu,0,0,rows,k);
461 BlockType A_2(lu,0,k+bs,rows,tsize);
462 BlockType A11(lu,k,k,bs,bs);
463 BlockType A12(lu,k,k+bs,bs,tsize);
464 BlockType A21(lu,k+bs,k,trows,bs);
465 BlockType A22(lu,k+bs,k+bs,trows,tsize);
467 PivIndex nb_transpositions_in_panel;
470 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
471 row_transpositions+k, nb_transpositions_in_panel, 16);
472 if(ret>=0 && first_zero_pivot==-1)
473 first_zero_pivot = k+ret;
475 nb_transpositions += nb_transpositions_in_panel;
477 for(
Index i=k; i<k+bs; ++i)
479 Index piv = (row_transpositions[i] += internal::convert_index<PivIndex>(k));
480 A_0.row(i).swap(A_0.row(piv));
486 for(
Index i=k;i<k+bs; ++i)
487 A_2.row(i).swap(A_2.row(row_transpositions[i]));
490 A11.template triangularView<UnitLower>().solveInPlace(A12);
492 A22.noalias() -= A21 * A12;
495 return first_zero_pivot;
501 template<
typename MatrixType,
typename TranspositionType>
502 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions,
typename TranspositionType::StorageIndex& nb_transpositions)
504 eigen_assert(lu.cols() == row_transpositions.size());
505 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
509 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
514 template<
typename MatrixType>
517 check_template_parameters();
522 m_l1_norm = m_lu.cwiseAbs().colwise().sum().maxCoeff();
524 eigen_assert(m_lu.rows() == m_lu.cols() &&
"PartialPivLU is only for square (and moreover invertible) matrices");
527 m_rowsTranspositions.resize(size);
529 typename TranspositionType::StorageIndex nb_transpositions;
530 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
531 m_det_p = (nb_transpositions%2) ? -1 : 1;
533 m_p = m_rowsTranspositions;
535 m_isInitialized =
true;
538 template<
typename MatrixType>
541 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
542 return Scalar(m_det_p) * m_lu.diagonal().prod();
548 template<
typename MatrixType>
551 eigen_assert(m_isInitialized &&
"LU is not initialized.");
553 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
554 * m_lu.template triangularView<Upper>();
557 res = m_p.inverse() * res;
567 template<
typename DstXprType,
typename MatrixType>
568 struct Assignment<DstXprType, Inverse<PartialPivLU<MatrixType> >, internal::assign_op<typename DstXprType::Scalar,typename PartialPivLU<MatrixType>::Scalar>, Dense2Dense>
572 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<typename DstXprType::Scalar,typename LuType::Scalar> &)
574 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
587 template<
typename Derived>
602 template<
typename Derived>
611 #endif // EIGEN_PARTIALLU_H Definition: Constants.h:320
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:292
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:141
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:88
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:175
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:549
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:250
Namespace containing all symbols from the Eigen library.
Definition: Core:287
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
Definition: NumTraits.h:150
Derived & derived()
Definition: EigenBase.h:45
Eigen::Index Index
The interface type of indices.
Definition: EigenBase.h:38
const unsigned int RowMajorBit
Definition: Constants.h:61
Definition: EigenBase.h:29
Expression of the inverse of another expression.
Definition: Inverse.h:43
Index size() const
Definition: EigenBase.h:66
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:33
Scalar determinant() const
Definition: PartialPivLU.h:539
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:197
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:589
Definition: Eigen_Colamd.h:50
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:103
const PermutationType & permutationP() const
Definition: PartialPivLU.h:149
Definition: Constants.h:322
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:604
const int Dynamic
Definition: Constants.h:21
Pseudo expression representing a solving operation.
Definition: Solve.h:62
RealScalar rcond() const
Definition: PartialPivLU.h:184
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48