11 #ifndef EIGEN_PARTIALLU_H
12 #define EIGEN_PARTIALLU_H
47 template<
typename _MatrixType>
class PartialPivLU
51 typedef _MatrixType MatrixType;
53 RowsAtCompileTime = MatrixType::RowsAtCompileTime,
54 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
55 Options = MatrixType::Options,
56 MaxRowsAtCompileTime = MatrixType::MaxRowsAtCompileTime,
57 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
59 typedef typename MatrixType::Scalar Scalar;
60 typedef typename NumTraits<typename MatrixType::Scalar>::Real RealScalar;
61 typedef typename internal::traits<MatrixType>::StorageKind StorageKind;
62 typedef typename MatrixType::Index Index;
63 typedef PermutationMatrix<RowsAtCompileTime, MaxRowsAtCompileTime> PermutationType;
64 typedef Transpositions<RowsAtCompileTime, MaxRowsAtCompileTime> TranspositionType;
102 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
110 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
131 template<
typename Rhs>
132 inline const internal::solve_retval<PartialPivLU, Rhs>
135 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
136 return internal::solve_retval<PartialPivLU, Rhs>(*
this, b.derived());
146 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
inverse()
const
148 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
149 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
150 (*
this, MatrixType::Identity(m_lu.rows(), m_lu.cols()));
166 typename internal::traits<MatrixType>::Scalar
determinant()
const;
170 inline Index rows()
const {
return m_lu.rows(); }
171 inline Index cols()
const {
return m_lu.cols(); }
176 TranspositionType m_rowsTranspositions;
178 bool m_isInitialized;
181 template<
typename MatrixType>
185 m_rowsTranspositions(),
187 m_isInitialized(false)
191 template<
typename MatrixType>
195 m_rowsTranspositions(size),
197 m_isInitialized(false)
201 template<
typename MatrixType>
203 : m_lu(matrix.rows(), matrix.rows()),
205 m_rowsTranspositions(matrix.rows()),
207 m_isInitialized(false)
215 template<
typename Scalar,
int StorageOrder,
typename PivIndex>
216 struct partial_lu_impl
226 typedef typename MatrixType::RealScalar RealScalar;
227 typedef typename MatrixType::Index Index;
239 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
241 const Index rows = lu.rows();
242 const Index cols = lu.cols();
243 const Index size = (std::min)(rows,cols);
244 nb_transpositions = 0;
245 Index first_zero_pivot = -1;
246 for(Index k = 0; k < size; ++k)
248 Index rrows = rows-k-1;
249 Index rcols = cols-k-1;
251 Index row_of_biggest_in_col;
252 RealScalar biggest_in_corner
253 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
254 row_of_biggest_in_col += k;
256 row_transpositions[k] = PivIndex(row_of_biggest_in_col);
258 if(biggest_in_corner != RealScalar(0))
260 if(k != row_of_biggest_in_col)
262 lu.row(k).swap(lu.row(row_of_biggest_in_col));
268 lu.col(k).tail(rrows) /= lu.coeff(k,k);
270 else if(first_zero_pivot==-1)
274 first_zero_pivot = k;
278 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
280 return first_zero_pivot;
298 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
300 MapLU lu1(lu_data,StorageOrder==
RowMajor?rows:luStride,StorageOrder==
RowMajor?luStride:cols);
301 MatrixType lu(lu1,0,0,rows,cols);
303 const Index size = (std::min)(rows,cols);
308 return unblocked_lu(lu, row_transpositions, nb_transpositions);
316 blockSize = (blockSize/16)*16;
317 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
320 nb_transpositions = 0;
321 Index first_zero_pivot = -1;
322 for(Index k = 0; k < size; k+=blockSize)
324 Index bs = (std::min)(size-k,blockSize);
325 Index trows = rows - k - bs;
326 Index tsize = size - k - bs;
332 BlockType A_0(lu,0,0,rows,k);
333 BlockType A_2(lu,0,k+bs,rows,tsize);
334 BlockType A11(lu,k,k,bs,bs);
335 BlockType A12(lu,k,k+bs,bs,tsize);
336 BlockType A21(lu,k+bs,k,trows,bs);
337 BlockType A22(lu,k+bs,k+bs,trows,tsize);
339 PivIndex nb_transpositions_in_panel;
342 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
343 row_transpositions+k, nb_transpositions_in_panel, 16);
344 if(ret>=0 && first_zero_pivot==-1)
345 first_zero_pivot = k+ret;
347 nb_transpositions += nb_transpositions_in_panel;
349 for(Index i=k; i<k+bs; ++i)
351 Index piv = (row_transpositions[i] += k);
352 A_0.row(i).swap(A_0.row(piv));
358 for(Index i=k;i<k+bs; ++i)
359 A_2.row(i).swap(A_2.row(row_transpositions[i]));
362 A11.template triangularView<UnitLower>().solveInPlace(A12);
364 A22.noalias() -= A21 * A12;
367 return first_zero_pivot;
373 template<
typename MatrixType,
typename TranspositionType>
374 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions,
typename TranspositionType::Index& nb_transpositions)
376 eigen_assert(lu.cols() == row_transpositions.size());
377 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
381 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
386 template<
typename MatrixType>
387 PartialPivLU<MatrixType>& PartialPivLU<MatrixType>::compute(
const MatrixType& matrix)
390 eigen_assert(matrix.rows()<NumTraits<int>::highest());
394 eigen_assert(matrix.rows() == matrix.cols() &&
"PartialPivLU is only for square (and moreover invertible) matrices");
395 const Index size = matrix.rows();
397 m_rowsTranspositions.resize(size);
399 typename TranspositionType::Index nb_transpositions;
400 internal::partial_lu_inplace(m_lu, m_rowsTranspositions, nb_transpositions);
401 m_det_p = (nb_transpositions%2) ? -1 : 1;
403 m_p = m_rowsTranspositions;
405 m_isInitialized =
true;
409 template<
typename MatrixType>
412 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
413 return Scalar(m_det_p) * m_lu.diagonal().prod();
419 template<
typename MatrixType>
422 eigen_assert(m_isInitialized &&
"LU is not initialized.");
424 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
425 * m_lu.template triangularView<Upper>();
428 res = m_p.inverse() * res;
437 template<
typename _MatrixType,
typename Rhs>
439 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
443 template<typename Dest>
void evalTo(Dest& dst)
const
452 eigen_assert(rhs().rows() == dec().matrixLU().rows());
455 dst = dec().permutationP() * rhs();
458 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
461 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
475 template<
typename Derived>
476 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
482 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
491 template<
typename Derived>
501 #endif // EIGEN_PARTIALLU_H
const PartialPivLU< PlainObject > lu() const
Definition: PartialPivLU.h:493
Definition: Constants.h:266
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:104
MatrixType reconstructedMatrix() const
Definition: PartialPivLU.h:420
const internal::solve_retval< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Definition: PartialPivLU.h:133
const PartialPivLU< PlainObject > partialPivLu() const
Definition: PartialPivLU.h:477
Definition: Constants.h:264
const MatrixType & matrixLU() const
Definition: PartialPivLU.h:100
LU decomposition of a matrix with partial pivoting, and related features.
Definition: ForwardDeclarations.h:217
internal::traits< MatrixType >::Scalar determinant() const
Definition: PartialPivLU.h:410
const internal::solve_retval< PartialPivLU, typename MatrixType::IdentityReturnType > inverse() const
Definition: PartialPivLU.h:146
Definition: Eigen_Colamd.h:54
Expression of a fixed-size or dynamic-size block.
Definition: Block.h:102
PartialPivLU()
Default Constructor.
Definition: PartialPivLU.h:182
const unsigned int RowMajorBit
Definition: Constants.h:53
const PermutationType & permutationP() const
Definition: PartialPivLU.h:108
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48