11 #ifndef EIGEN_PARTIALLU_H 12 #define EIGEN_PARTIALLU_H 47 template<
typename _MatrixType>
class PartialPivLU
131 template<
typename Rhs>
150 (*
this, MatrixType::Identity(
m_lu.rows(),
m_lu.cols()));
181 template<
typename MatrixType>
191 template<
typename MatrixType>
201 template<
typename MatrixType>
215 template<
typename Scalar,
int StorageOrder,
typename PivIndex>
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);
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>
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>
394 eigen_assert(matrix.rows() == matrix.cols() &&
"PartialPivLU is only for square (and moreover invertible) matrices");
395 const Index size = matrix.rows();
401 m_det_p = (nb_transpositions%2) ? -1 : 1;
409 template<
typename MatrixType>
419 template<
typename MatrixType>
424 MatrixType res =
m_lu.template triangularView<UnitLower>().toDenseMatrix()
425 *
m_lu.template triangularView<Upper>();
437 template<
typename _MatrixType,
typename Rhs>
443 template<typename Dest>
void evalTo(Dest& dst)
const 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>
482 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS 491 template<
typename Derived>
501 #endif // EIGEN_PARTIALLU_H
internal::traits< Derived >::Scalar Scalar
MatrixType::Scalar Scalar
PartialPivLU()
Default Constructor.
Block< MatrixType, Dynamic, Dynamic > BlockType
Map< Matrix< Scalar, Dynamic, Dynamic, StorageOrder > > MapLU
A matrix or vector expression mapping an existing array of data.
internal::traits< Derived >::Index Index
PartialPivLU & compute(const MatrixType &matrix)
LU decomposition of a matrix with partial pivoting, and related features.
iterative scaling algorithm to equilibrate rows and column norms in matrices
void partial_lu_inplace(MatrixType &lu, TranspositionType &row_transpositions, typename TranspositionType::Index &nb_transpositions)
Holds information about the various numeric (i.e. scalar) types allowed by Eigen. ...
const internal::solve_retval< PartialPivLU, Rhs > solve(const MatrixBase< Rhs > &b) const
Transpositions< RowsAtCompileTime, MaxRowsAtCompileTime > TranspositionType
MatrixType::RealScalar RealScalar
const unsigned int RowMajorBit
const PartialPivLU< PlainObject > partialPivLu() const
const PermutationType & permutationP() const
Provides a generic way to set and pass user-specified options.
const internal::solve_retval< PartialPivLU, typename MatrixType::IdentityReturnType > inverse() const
Transpose< PermutationBase > inverse() const
void rhs(const real_t *x, real_t *f)
PermutationMatrix< RowsAtCompileTime, MaxRowsAtCompileTime > PermutationType
Expression of a fixed-size or dynamic-size block.
static Index unblocked_lu(MatrixType &lu, PivIndex *row_transpositions, PivIndex &nb_transpositions)
Block< MapLU, Dynamic, Dynamic > MatrixType
internal::traits< MatrixType >::Scalar determinant() const
TranspositionType m_rowsTranspositions
NumTraits< typename MatrixType::Scalar >::Real RealScalar
#define EIGEN_MAKE_SOLVE_HELPERS(DecompositionType, Rhs)
MatrixType reconstructedMatrix() const
const MatrixType & matrixLU() const
internal::traits< MatrixType >::StorageKind StorageKind
static Index blocked_lu(Index rows, Index cols, Scalar *lu_data, Index luStride, PivIndex *row_transpositions, PivIndex &nb_transpositions, Index maxBlockSize=256)
Base class for all dense matrices, vectors, and expressions.
IndicesType::Scalar Index