12 #ifndef EIGEN_SPARSE_LU_H
13 #define EIGEN_SPARSE_LU_H
17 template <
typename _MatrixType,
typename _OrderingType = COLAMDOrdering<
typename _MatrixType::StorageIndex> >
class SparseLU;
18 template <
typename MappedSparseMatrixType>
struct SparseLUMatrixLReturnType;
19 template <
typename MatrixLType,
typename MatrixUType>
struct SparseLUMatrixUReturnType;
72 template <
typename _MatrixType,
typename _OrderingType>
73 class SparseLU :
public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex>
77 using APIBase::m_isInitialized;
79 using APIBase::_solve_impl;
81 typedef _MatrixType MatrixType;
82 typedef _OrderingType OrderingType;
83 typedef typename MatrixType::Scalar Scalar;
84 typedef typename MatrixType::RealScalar RealScalar;
85 typedef typename MatrixType::StorageIndex StorageIndex;
87 typedef internal::MappedSuperNodalMatrix<Scalar, StorageIndex> SCMatrix;
91 typedef internal::SparseLUImpl<Scalar, StorageIndex> Base;
94 SparseLU():m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
98 explicit SparseLU(
const MatrixType& matrix):m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
110 void factorize (
const MatrixType& matrix);
111 void simplicialfactorize(
const MatrixType& matrix);
125 inline Index rows()
const {
return m_mat.
rows(); }
126 inline Index cols()
const {
return m_mat.
cols(); }
130 m_symmetricmode = sym;
139 SparseLUMatrixLReturnType<SCMatrix>
matrixL()
const
141 return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
149 SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >
matrixU()
const
151 return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >(m_Lstore, m_Ustore);
173 m_diagpivotthresh = thresh;
176 #ifdef EIGEN_PARSED_BY_DOXYGEN
183 template<
typename Rhs>
185 #endif // EIGEN_PARSED_BY_DOXYGEN
197 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
209 template<
typename Rhs,
typename Dest>
212 Dest& X(X_base.derived());
213 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first");
215 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
219 X.resize(B.rows(),B.cols());
222 for(Index j = 0; j < B.cols(); ++j)
226 this->
matrixL().solveInPlace(X);
227 this->
matrixU().solveInPlace(X);
230 for (Index j = 0; j < B.cols(); ++j)
249 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
251 Scalar det = Scalar(1.);
254 for (Index j = 0; j < this->cols(); ++j)
256 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
260 det *= abs(it.value());
281 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
282 Scalar det = Scalar(0.);
283 for (Index j = 0; j < this->cols(); ++j)
285 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
287 if(it.row() < j)
continue;
290 det += log(abs(it.value()));
304 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
309 for (Index j = 0; j < this->cols(); ++j)
311 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
317 else if(it.value()==0)
323 return det * m_detPermR * m_detPermC;
332 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
334 Scalar det = Scalar(1.);
337 for (Index j = 0; j < this->cols(); ++j)
339 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
348 return (m_detPermR * m_detPermC) > 0 ? det : -det;
353 void initperfvalues()
355 m_perfv.panel_size = 16;
357 m_perfv.maxsuper = 128;
360 m_perfv.fillfactor = 20;
365 bool m_factorizationIsOk;
367 std::string m_lastError;
370 MappedSparseMatrix<Scalar,ColMajor,StorageIndex> m_Ustore;
371 PermutationType m_perm_c;
372 PermutationType m_perm_r ;
375 typename Base::GlobalLU_t m_glu;
378 bool m_symmetricmode;
380 internal::perfvalues m_perfv;
381 RealScalar m_diagpivotthresh;
382 Index m_nnzL, m_nnzU;
383 Index m_detPermR, m_detPermC;
386 SparseLU (
const SparseLU& );
403 template <
typename MatrixType,
typename OrderingType>
421 ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?
const_cast<StorageIndex*
>(mat.outerIndexPtr()):0);
424 if(!mat.isCompressed())
425 IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
428 for (Index i = 0; i < mat.cols(); i++)
430 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
431 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
436 IndexVector firstRowElt;
437 internal::coletree(m_mat, m_etree,firstRowElt);
440 if (!m_symmetricmode) {
441 IndexVector post, iwork;
443 internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post);
447 Index m = m_mat.cols();
449 for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
453 PermutationType post_perm(m);
454 for (Index i = 0; i < m; i++)
455 post_perm.
indices()(i) = post(i);
458 if(m_perm_c.size()) {
459 m_perm_c = post_perm * m_perm_c;
464 m_analysisIsOk =
true;
488 template <
typename MatrixType,
typename OrderingType>
491 using internal::emptyIdxLU;
492 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
493 eigen_assert((matrix.rows() == matrix.cols()) &&
"Only for squared matrices");
495 typedef typename IndexVector::Scalar StorageIndex;
497 m_isInitialized =
true;
507 const StorageIndex * outerIndexPtr;
508 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
511 StorageIndex* outerIndexPtr_t =
new StorageIndex[matrix.cols()+1];
512 for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
513 outerIndexPtr = outerIndexPtr_t;
515 for (Index i = 0; i < matrix.cols(); i++)
517 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
518 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
520 if(!matrix.isCompressed())
delete[] outerIndexPtr;
524 m_perm_c.resize(matrix.cols());
525 for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
528 Index m = m_mat.rows();
529 Index n = m_mat.cols();
530 Index nnz = m_mat.nonZeros();
531 Index maxpanel = m_perfv.panel_size * m;
534 Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
537 m_lastError =
"UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
538 m_factorizationIsOk =
false;
543 IndexVector segrep(m); segrep.
setZero();
544 IndexVector parent(m); parent.
setZero();
545 IndexVector xplore(m); xplore.
setZero();
546 IndexVector repfnz(maxpanel);
547 IndexVector panel_lsub(maxpanel);
548 IndexVector xprune(n); xprune.
setZero();
549 IndexVector marker(m*internal::LUNoMarker); marker.
setZero();
558 tempv.
setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, m) );
561 PermutationType iperm_c(m_perm_c.inverse());
564 IndexVector relax_end(n);
565 if ( m_symmetricmode ==
true )
566 Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
568 Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
572 m_perm_r.indices().setConstant(-1);
576 m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
577 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
583 IndexVector panel_histo(n);
589 for (jcol = 0; jcol < n; )
592 Index panel_size = m_perfv.panel_size;
593 for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
595 if (relax_end(k) != emptyIdxLU)
597 panel_size = k - jcol;
602 panel_size = n - jcol;
605 Base::panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu);
608 Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
611 for ( jj = jcol; jj< jcol + panel_size; jj++)
619 info = Base::column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu);
622 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
624 m_factorizationIsOk =
false;
630 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
633 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
635 m_factorizationIsOk =
false;
640 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
643 m_lastError =
"UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
645 m_factorizationIsOk =
false;
650 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
653 m_lastError =
"THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
654 std::ostringstream returnInfo;
656 m_lastError += returnInfo.str();
658 m_factorizationIsOk =
false;
664 if (pivrow != jj) m_detPermR = -m_detPermR;
667 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
670 for (i = 0; i < nseg; i++)
673 repfnz_k(irep) = emptyIdxLU;
679 m_detPermR = m_perm_r.determinant();
680 m_detPermC = m_perm_c.determinant();
683 Base::countnz(n, m_nnzL, m_nnzU, m_glu);
685 Base::fixupL(n, m_perm_r.indices(), m_glu);
688 m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
693 m_factorizationIsOk =
true;
696 template<
typename MappedSupernodalType>
697 struct SparseLUMatrixLReturnType : internal::no_assignment_operator
699 typedef typename MappedSupernodalType::Scalar Scalar;
700 explicit SparseLUMatrixLReturnType(
const MappedSupernodalType& mapL) : m_mapL(mapL)
702 Index rows() {
return m_mapL.rows(); }
703 Index cols() {
return m_mapL.cols(); }
704 template<
typename Dest>
705 void solveInPlace( MatrixBase<Dest> &X)
const
707 m_mapL.solveInPlace(X);
709 const MappedSupernodalType& m_mapL;
712 template<
typename MatrixLType,
typename MatrixUType>
713 struct SparseLUMatrixUReturnType : internal::no_assignment_operator
715 typedef typename MatrixLType::Scalar Scalar;
716 explicit SparseLUMatrixUReturnType(
const MatrixLType& mapL,
const MatrixUType& mapU)
717 : m_mapL(mapL),m_mapU(mapU)
719 Index rows() {
return m_mapL.rows(); }
720 Index cols() {
return m_mapL.cols(); }
722 template<
typename Dest>
void solveInPlace(MatrixBase<Dest> &X)
const
724 Index nrhs = X.cols();
727 for (Index k = m_mapL.nsuper(); k >= 0; k--)
729 Index fsupc = m_mapL.supToCol()[k];
730 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc];
731 Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
732 Index luptr = m_mapL.colIndexPtr()[fsupc];
736 for (Index j = 0; j < nrhs; j++)
738 X(fsupc, j) /= m_mapL.valuePtr()[luptr];
743 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
744 Map< Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
745 U = A.template triangularView<Upper>().solve(U);
748 for (Index j = 0; j < nrhs; ++j)
750 for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
752 typename MatrixUType::InnerIterator it(m_mapU, jcol);
755 Index irow = it.index();
756 X(irow, j) -= X(jcol, j) * it.value();
762 const MatrixLType& m_mapL;
763 const MatrixUType& m_mapU;
ColXpr col(Index i)
Definition: DenseBase.h:778
Index cols() const
Definition: SparseMatrix.h:132
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:158
A base class for sparse solvers.
Definition: SparseSolverBase.h:53
std::string lastErrorMessage() const
Definition: SparseLU.h:204
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:489
const unsigned int RowMajorBit
Definition: Constants.h:53
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:252
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:195
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:171
Scalar absDeterminant()
Definition: SparseLU.h:246
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, StorageIndex > > matrixU() const
Definition: SparseLU.h:149
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:520
Expression of a fixed-size or dynamic-size sub-vector.
Definition: ForwardDeclarations.h:87
void compute(const MatrixType &matrix)
Definition: SparseLU.h:117
const PermutationType & colsPermutation() const
Definition: SparseLU.h:166
void isSymmetric(bool sym)
Definition: SparseLU.h:128
Definition: Constants.h:426
const Solve< SparseLU, Rhs > solve(const MatrixBase< Rhs > &B) const
Scalar determinant()
Definition: SparseLU.h:330
Definition: Constants.h:424
const IndicesType & indices() const
Definition: PermutationMatrix.h:390
Scalar logAbsDeterminant() const
Definition: SparseLU.h:276
Pseudo expression representing a solving operation.
Definition: Solve.h:63
TransposeReturnType inverse() const
Definition: PermutationMatrix.h:198
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:404
ComputationInfo
Definition: Constants.h:422
Scalar signDeterminant()
Definition: SparseLU.h:302
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:139
Index rows() const
Definition: SparseMatrix.h:130
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:352