11 #ifndef EIGEN_PARTIALLU_H
12 #define EIGEN_PARTIALLU_H
17 template<
typename _MatrixType>
struct traits<PartialPivLU<_MatrixType> >
20 typedef traits<_MatrixType> BaseTraits;
23 CoeffReadCost = Dynamic
60 template<
typename _MatrixType>
class PartialPivLU
64 typedef _MatrixType MatrixType;
66 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
67 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
68 Options = MatrixType::Options,
69 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
70 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
72 typedef typename MatrixType::Scalar Scalar;
73 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
74 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
76 typedef typename MatrixType::StorageIndex StorageIndex;
77 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
78 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
79 typedef typename MatrixType::PlainObject PlainObject;
117 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
125 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
146 template<
typename Rhs>
150 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
163 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
180 typename internal::traits<MatrixType>::Scalar
determinant()
const;
184 inline Index rows()
const {
return m_lu.rows(); }
185 inline Index cols()
const {
return m_lu.cols(); }
187 #ifndef EIGEN_PARSED_BY_DOXYGEN
188 template<
typename RhsType,
typename DstType>
190 void _solve_impl(
const RhsType &rhs, DstType &dst)
const {
198 eigen_assert(rhs.rows() == m_lu.rows());
204 m_lu.template triangularView<UnitLower>().solveInPlace(dst);
207 m_lu.template triangularView<Upper>().solveInPlace(dst);
213 static void check_template_parameters()
215 EIGEN_STATIC_ASSERT_NON_INTEGER(Scalar);
220 TranspositionType m_rowsTranspositions;
222 bool m_isInitialized;
225 template<
typename MatrixType>
229 m_rowsTranspositions(),
231 m_isInitialized(false)
235 template<
typename MatrixType>
239 m_rowsTranspositions(size),
241 m_isInitialized(false)
245 template<
typename MatrixType>
247 : m_lu(matrix.rows(), matrix.rows()),
249 m_rowsTranspositions(matrix.rows()),
251 m_isInitialized(false)
259 template<
typename Scalar,
int StorageOrder,
typename PivIndex>
260 struct partial_lu_impl
270 typedef typename MatrixType::RealScalar RealScalar;
282 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
284 typedef scalar_score_coeff_op<Scalar> Scoring;
285 typedef typename Scoring::result_type Score;
286 const Index rows = lu.rows();
287 const Index cols = lu.cols();
288 const Index size = (std::min)(rows,cols);
289 nb_transpositions = 0;
290 Index first_zero_pivot = -1;
291 for(Index k = 0; k < size; ++k)
293 Index rrows = rows-k-1;
294 Index rcols = cols-k-1;
296 Index row_of_biggest_in_col;
297 Score biggest_in_corner
298 = lu.col(k).tail(rows-k).unaryExpr(Scoring()).maxCoeff(&row_of_biggest_in_col);
299 row_of_biggest_in_col += k;
301 row_transpositions[k] = PivIndex(row_of_biggest_in_col);
303 if(biggest_in_corner != Score(0))
305 if(k != row_of_biggest_in_col)
307 lu.row(k).swap(lu.row(row_of_biggest_in_col));
313 lu.col(k).tail(rrows) /= lu.coeff(k,k);
315 else if(first_zero_pivot==-1)
319 first_zero_pivot = k;
323 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
325 return first_zero_pivot;
343 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
345 MapLU lu1(lu_data,StorageOrder==
RowMajor?rows:luStride,StorageOrder==
RowMajor?luStride:cols);
346 MatrixType lu(lu1,0,0,rows,cols);
348 const Index size = (std::min)(rows,cols);
353 return unblocked_lu(lu, row_transpositions, nb_transpositions);
361 blockSize = (blockSize/16)*16;
362 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
365 nb_transpositions = 0;
366 Index first_zero_pivot = -1;
367 for(Index k = 0; k < size; k+=blockSize)
369 Index bs = (std::min)(size-k,blockSize);
370 Index trows = rows - k - bs;
371 Index tsize = size - k - bs;
377 BlockType A_0(lu,0,0,rows,k);
378 BlockType A_2(lu,0,k+bs,rows,tsize);
379 BlockType A11(lu,k,k,bs,bs);
380 BlockType A12(lu,k,k+bs,bs,tsize);
381 BlockType A21(lu,k+bs,k,trows,bs);
382 BlockType A22(lu,k+bs,k+bs,trows,tsize);
384 PivIndex nb_transpositions_in_panel;
387 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
388 row_transpositions+k, nb_transpositions_in_panel, 16);
389 if(ret>=0 && first_zero_pivot==-1)
390 first_zero_pivot = k+ret;
392 nb_transpositions += nb_transpositions_in_panel;
394 for(Index i=k; i<k+bs; ++i)
396 Index piv = (row_transpositions[i] += k);
397 A_0.row(i).swap(A_0.row(piv));
403 for(Index i=k;i<k+bs; ++i)
404 A_2.row(i).swap(A_2.row(row_transpositions[i]));
407 A11.template triangularView<UnitLower>().solveInPlace(A12);
409 A22.noalias() -= A21 * A12;
412 return first_zero_pivot;
418 template<
typename MatrixType,
typename TranspositionType>
419 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions,
typename TranspositionType::StorageIndex& nb_transpositions)
421 eigen_assert(lu.cols() == row_transpositions.size());
422 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
426 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
431 template<
typename MatrixType>
432 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(
const MatrixType& matrix)
434 check_template_parameters();
437 eigen_assert(matrix.rows()<NumTraits<int>::highest());
441 eigen_assert(matrix.rows() == matrix.cols() &&
"PartialPivLU is only for square (and moreover invertible) matrices");
442 const Index size = matrix.rows();
444 m_rowsTranspositions.resize(size);
446 typename TranspositionType::StorageIndex nb_transpositions;
447 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
448 m_det_p = (nb_transpositions%2) ? -1 : 1;
450 m_p = m_rowsTranspositions;
452 m_isInitialized =
true;
456 template<
typename MatrixType>
459 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
460 return Scalar(m_det_p) * m_lu.diagonal().prod();
466 template<
typename MatrixType>
469 eigen_assert(m_isInitialized &&
"LU is not initialized.");
471 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
472 * m_lu.template triangularView<Upper>();
475 res = m_p.inverse() * res;
485 template<
typename DstXprType,
typename MatrixType,
typename Scalar>
490 static void run(DstXprType &dst,
const SrcXprType &src,
const internal::assign_op<Scalar> &)
492 dst = src.nestedExpression().solve(MatrixType::Identity(src.rows(), src.cols()));
506 template<
typename Derived>
507 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
523 template<
typename Derived>
533 #endif // EIGEN_PARTIALLU_H
Definition: Constants.h:314
const Inverse< PartialPivLU > inverse() const
Definition: PartialPivLU.h:161
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:226
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:89
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:247
const unsigned int RowMajorBit
Definition: Constants.h:53
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:508
Expression of the inverse of another expression.
Definition: Inverse.h:45
const PermutationType & permutationP() const
Definition: PartialPivLU.h:123
Definition: Eigen_Colamd.h:54
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:104
internal::traits< MatrixType >::Scalar determinant() const
Definition: PartialPivLU.h:457
Definition: Constants.h:312
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:467
Pseudo expression representing a solving operation.
Definition: Solve.h:63
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:115
const Solve< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:148
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:525