26 #ifndef EIGEN_PARTIALLU_H
27 #define EIGEN_PARTIALLU_H
74 typedef typename MatrixType::Scalar
Scalar;
76 typedef typename internal::traits<MatrixType>::StorageKind
StorageKind;
77 typedef typename MatrixType::Index
Index;
146 template<
typename Rhs>
147 inline const internal::solve_retval<PartialPivLU, Rhs>
151 return internal::solve_retval<PartialPivLU, Rhs>(*
this, b.derived());
161 inline const internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
inverse()
const
164 return internal::solve_retval<PartialPivLU,typename MatrixType::IdentityReturnType>
165 (*
this, MatrixType::Identity(
m_lu.rows(),
m_lu.cols()));
181 typename internal::traits<MatrixType>::Scalar
determinant()
const;
196 template<
typename MatrixType>
200 m_rowsTranspositions(),
202 m_isInitialized(false)
206 template<
typename MatrixType>
210 m_rowsTranspositions(size),
212 m_isInitialized(false)
216 template<
typename MatrixType>
218 : m_lu(matrix.rows(), matrix.rows()),
220 m_rowsTranspositions(matrix.rows()),
222 m_isInitialized(false)
230 template<
typename Scalar,
int StorageOrder,
typename PivIndex>
231 struct partial_lu_impl
241 typedef typename MatrixType::RealScalar RealScalar;
242 typedef typename MatrixType::Index Index;
254 static Index unblocked_lu(MatrixType& lu, PivIndex* row_transpositions, PivIndex& nb_transpositions)
256 const Index rows = lu.
rows();
257 const Index cols = lu.cols();
258 const Index size = (std::min)(rows,cols);
259 nb_transpositions = 0;
260 int first_zero_pivot = -1;
261 for(Index k = 0; k < size; ++k)
263 Index rrows = rows-k-1;
264 Index rcols = cols-k-1;
266 Index row_of_biggest_in_col;
267 RealScalar biggest_in_corner
268 = lu.col(k).tail(rows-k).cwiseAbs().maxCoeff(&row_of_biggest_in_col);
269 row_of_biggest_in_col += k;
271 row_transpositions[k] = row_of_biggest_in_col;
273 if(biggest_in_corner != RealScalar(0))
275 if(k != row_of_biggest_in_col)
277 lu.row(k).swap(lu.row(row_of_biggest_in_col));
283 lu.col(k).tail(rrows) /= lu.coeff(k,k);
285 else if(first_zero_pivot==-1)
289 first_zero_pivot = k;
293 lu.bottomRightCorner(rrows,rcols).noalias() -= lu.col(k).tail(rrows) * lu.row(k).tail(rcols);
295 return first_zero_pivot;
313 static Index blocked_lu(Index rows, Index cols, Scalar* lu_data, Index luStride, PivIndex* row_transpositions, PivIndex& nb_transpositions, Index maxBlockSize=256)
315 MapLU lu1(lu_data,StorageOrder==
RowMajor?rows:luStride,StorageOrder==
RowMajor?luStride:cols);
316 MatrixType lu(lu1,0,0,rows,cols);
318 const Index size = (std::min)(rows,cols);
323 return unblocked_lu(lu, row_transpositions, nb_transpositions);
331 blockSize = (blockSize/16)*16;
332 blockSize = (std::min)((std::max)(blockSize,Index(8)), maxBlockSize);
335 nb_transpositions = 0;
336 int first_zero_pivot = -1;
337 for(Index k = 0; k < size; k+=blockSize)
339 Index bs = (std::min)(size-k,blockSize);
340 Index trows = rows - k - bs;
341 Index tsize = size - k - bs;
347 BlockType A_0(lu,0,0,rows,k);
348 BlockType A_2(lu,0,k+bs,rows,tsize);
349 BlockType A11(lu,k,k,bs,bs);
350 BlockType A12(lu,k,k+bs,bs,tsize);
351 BlockType A21(lu,k+bs,k,trows,bs);
352 BlockType A22(lu,k+bs,k+bs,trows,tsize);
354 PivIndex nb_transpositions_in_panel;
357 Index ret = blocked_lu(trows+bs, bs, &lu.coeffRef(k,k), luStride,
358 row_transpositions+k, nb_transpositions_in_panel, 16);
359 if(ret>=0 && first_zero_pivot==-1)
360 first_zero_pivot = k+ret;
362 nb_transpositions += nb_transpositions_in_panel;
364 for(Index i=k; i<k+bs; ++i)
366 Index piv = (row_transpositions[i] += k);
367 A_0.row(i).swap(A_0.row(piv));
373 for(Index i=k;i<k+bs; ++i)
374 A_2.row(i).swap(A_2.row(row_transpositions[i]));
377 A11.template triangularView<UnitLower>().solveInPlace(A12);
379 A22.noalias() -= A21 * A12;
382 return first_zero_pivot;
388 template<
typename MatrixType,
typename TranspositionType>
389 void partial_lu_inplace(MatrixType& lu, TranspositionType& row_transpositions,
typename TranspositionType::Index& nb_transpositions)
392 eigen_assert((&row_transpositions.coeffRef(1)-&row_transpositions.coeffRef(0)) == 1);
396 ::blocked_lu(lu.rows(), lu.cols(), &lu.coeffRef(0,0), lu.outerStride(), &row_transpositions.coeffRef(0), nb_transpositions);
401 template<
typename MatrixType>
406 eigen_assert(matrix.rows() == matrix.cols() &&
"PartialPivLU is only for square (and moreover invertible) matrices");
407 const Index size = matrix.rows();
409 m_rowsTranspositions.resize(size);
413 m_det_p = (nb_transpositions%2) ? -1 : 1;
415 m_p = m_rowsTranspositions;
417 m_isInitialized =
true;
421 template<
typename MatrixType>
424 eigen_assert(m_isInitialized &&
"PartialPivLU is not initialized.");
425 return Scalar(m_det_p) * m_lu.diagonal().prod();
431 template<
typename MatrixType>
434 eigen_assert(m_isInitialized &&
"LU is not initialized.");
436 MatrixType res = m_lu.template triangularView<UnitLower>().toDenseMatrix()
437 * m_lu.template triangularView<Upper>();
440 res = m_p.inverse() * res;
449 template<
typename _MatrixType,
typename Rhs>
451 : solve_retval_base<PartialPivLU<_MatrixType>, Rhs>
455 template<typename Dest>
void evalTo(Dest& dst)
const
467 dst = dec().permutationP() * rhs();
470 dec().matrixLU().template triangularView<UnitLower>().solveInPlace(dst);
473 dec().matrixLU().template triangularView<Upper>().solveInPlace(dst);
487 template<
typename Derived>
488 inline const PartialPivLU<typename MatrixBase<Derived>::PlainObject>
494 #if EIGEN2_SUPPORT_STAGE > STAGE20_RESOLVE_API_CONFLICTS
503 template<
typename Derived>
513 #endif // EIGEN_PARTIALLU_H