12 #ifndef EIGEN_SPARSE_LU_H
13 #define EIGEN_SPARSE_LU_H
17 template <
typename _MatrixType,
typename _OrderingType = COLAMDOrdering<
typename _MatrixType::Index> >
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 internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index>
76 typedef _MatrixType MatrixType;
77 typedef _OrderingType OrderingType;
78 typedef typename MatrixType::Scalar Scalar;
79 typedef typename MatrixType::RealScalar RealScalar;
80 typedef typename MatrixType::Index Index;
82 typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix;
86 typedef internal::SparseLUImpl<Scalar, Index> Base;
89 SparseLU():m_isInitialized(true),m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
93 SparseLU(
const MatrixType& matrix):m_isInitialized(true),m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
105 void factorize (
const MatrixType& matrix);
106 void simplicialfactorize(
const MatrixType& matrix);
120 inline Index rows()
const {
return m_mat.
rows(); }
121 inline Index cols()
const {
return m_mat.
cols(); }
125 m_symmetricmode = sym;
134 SparseLUMatrixLReturnType<SCMatrix>
matrixL()
const
136 return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
144 SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> >
matrixU()
const
146 return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore);
168 m_diagpivotthresh = thresh;
177 template<
typename Rhs>
180 eigen_assert(m_factorizationIsOk &&
"SparseLU is not initialized.");
181 eigen_assert(rows()==B.rows()
182 &&
"SparseLU::solve(): invalid number of rows of the right hand side matrix B");
183 return internal::solve_retval<SparseLU, Rhs>(*
this, B.derived());
190 template<
typename Rhs>
193 eigen_assert(m_factorizationIsOk &&
"SparseLU is not initialized.");
194 eigen_assert(rows()==B.
rows()
195 &&
"SparseLU::solve(): invalid number of rows of the right hand side matrix B");
196 return internal::sparse_solve_retval<SparseLU, Rhs>(*
this, B.
derived());
209 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
221 template<
typename Rhs,
typename Dest>
224 Dest& X(X_base.derived());
225 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first");
227 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
231 X.resize(B.rows(),B.cols());
234 for(Index j = 0; j < B.cols(); ++j)
238 this->
matrixL().solveInPlace(X);
239 this->
matrixU().solveInPlace(X);
242 for (Index j = 0; j < B.cols(); ++j)
260 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
262 Scalar det = Scalar(1.);
265 for (Index j = 0; j < this->cols(); ++j)
267 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
269 if(it.row() < j)
continue;
272 det *= (std::abs)(it.value());
290 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
291 Scalar det = Scalar(0.);
292 for (Index j = 0; j < this->cols(); ++j)
294 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
296 if(it.row() < j)
continue;
299 det += (std::log)((std::abs)(it.value()));
313 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
314 return Scalar(m_detPermR);
319 void initperfvalues()
321 m_perfv.panel_size = 1;
323 m_perfv.maxsuper = 128;
326 m_perfv.fillfactor = 20;
331 bool m_isInitialized;
332 bool m_factorizationIsOk;
334 std::string m_lastError;
337 MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore;
338 PermutationType m_perm_c;
339 PermutationType m_perm_r ;
342 typename Base::GlobalLU_t m_glu;
345 bool m_symmetricmode;
347 internal::perfvalues<Index> m_perfv;
348 RealScalar m_diagpivotthresh;
349 Index m_nnzL, m_nnzU;
353 SparseLU (
const SparseLU& );
370 template <
typename MatrixType,
typename OrderingType>
382 if (m_perm_c.size()) {
385 const Index * outerIndexPtr;
386 if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr();
389 Index *outerIndexPtr_t =
new Index[mat.cols()+1];
390 for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
391 outerIndexPtr = outerIndexPtr_t;
393 for (Index i = 0; i < mat.cols(); i++)
395 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
396 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
398 if(!mat.isCompressed())
delete[] outerIndexPtr;
401 IndexVector firstRowElt;
405 if (!m_symmetricmode) {
406 IndexVector post, iwork;
412 Index m = m_mat.cols();
414 for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
418 PermutationType post_perm(m);
419 for (Index i = 0; i < m; i++)
420 post_perm.
indices()(i) = post(i);
423 if(m_perm_c.size()) {
424 m_perm_c = post_perm * m_perm_c;
429 m_analysisIsOk =
true;
453 template <
typename MatrixType,
typename OrderingType>
456 using internal::emptyIdxLU;
457 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
458 eigen_assert((matrix.rows() == matrix.cols()) &&
"Only for squared matrices");
460 typedef typename IndexVector::Scalar Index;
470 const Index * outerIndexPtr;
471 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
474 Index* outerIndexPtr_t =
new Index[matrix.cols()+1];
475 for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
476 outerIndexPtr = outerIndexPtr_t;
478 for (Index i = 0; i < matrix.cols(); i++)
480 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
481 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
483 if(!matrix.isCompressed())
delete[] outerIndexPtr;
487 m_perm_c.resize(matrix.cols());
488 for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
491 Index m = m_mat.rows();
492 Index n = m_mat.cols();
493 Index nnz = m_mat.nonZeros();
494 Index maxpanel = m_perfv.panel_size * m;
497 Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
500 m_lastError =
"UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
501 m_factorizationIsOk =
false;
506 IndexVector segrep(m); segrep.
setZero();
507 IndexVector parent(m); parent.
setZero();
508 IndexVector xplore(m); xplore.
setZero();
509 IndexVector repfnz(maxpanel);
510 IndexVector panel_lsub(maxpanel);
511 IndexVector xprune(n); xprune.
setZero();
512 IndexVector marker(m*internal::LUNoMarker); marker.
setZero();
521 tempv.
setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, m) );
524 PermutationType iperm_c(m_perm_c.inverse());
527 IndexVector relax_end(n);
528 if ( m_symmetricmode ==
true )
529 Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
531 Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
535 m_perm_r.indices().setConstant(-1);
539 m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
540 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
546 IndexVector panel_histo(n);
552 for (jcol = 0; jcol < n; )
555 Index panel_size = m_perfv.panel_size;
556 for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
558 if (relax_end(k) != emptyIdxLU)
560 panel_size = k - jcol;
565 panel_size = n - jcol;
568 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);
571 Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
574 for ( jj = jcol; jj< jcol + panel_size; jj++)
582 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);
585 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
587 m_factorizationIsOk =
false;
593 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
596 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
598 m_factorizationIsOk =
false;
603 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
606 m_lastError =
"UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
608 m_factorizationIsOk =
false;
613 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
616 m_lastError =
"THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
617 std::ostringstream returnInfo;
619 m_lastError += returnInfo.str();
621 m_factorizationIsOk =
false;
626 if (pivrow != jj) m_detPermR *= -1;
629 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
632 for (i = 0; i < nseg; i++)
635 repfnz_k(irep) = emptyIdxLU;
642 Base::countnz(n, m_nnzL, m_nnzU, m_glu);
644 Base::fixupL(n, m_perm_r.indices(), m_glu);
647 m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
652 m_factorizationIsOk =
true;
655 template<
typename MappedSupernodalType>
656 struct SparseLUMatrixLReturnType : internal::no_assignment_operator
658 typedef typename MappedSupernodalType::Index Index;
659 typedef typename MappedSupernodalType::Scalar Scalar;
660 SparseLUMatrixLReturnType(
const MappedSupernodalType& mapL) : m_mapL(mapL)
662 Index rows() {
return m_mapL.rows(); }
663 Index cols() {
return m_mapL.cols(); }
664 template<
typename Dest>
665 void solveInPlace( MatrixBase<Dest> &X)
const
667 m_mapL.solveInPlace(X);
669 const MappedSupernodalType& m_mapL;
672 template<
typename MatrixLType,
typename MatrixUType>
673 struct SparseLUMatrixUReturnType : internal::no_assignment_operator
675 typedef typename MatrixLType::Index Index;
676 typedef typename MatrixLType::Scalar Scalar;
677 SparseLUMatrixUReturnType(
const MatrixLType& mapL,
const MatrixUType& mapU)
678 : m_mapL(mapL),m_mapU(mapU)
680 Index rows() {
return m_mapL.rows(); }
681 Index cols() {
return m_mapL.cols(); }
683 template<
typename Dest>
void solveInPlace(MatrixBase<Dest> &X)
const
685 Index nrhs = X.cols();
688 for (Index k = m_mapL.nsuper(); k >= 0; k--)
690 Index fsupc = m_mapL.supToCol()[k];
691 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc];
692 Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
693 Index luptr = m_mapL.colIndexPtr()[fsupc];
697 for (Index j = 0; j < nrhs; j++)
699 X(fsupc, j) /= m_mapL.valuePtr()[luptr];
704 Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
705 Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
706 U = A.template triangularView<Upper>().solve(U);
709 for (Index j = 0; j < nrhs; ++j)
711 for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
713 typename MatrixUType::InnerIterator it(m_mapU, jcol);
716 Index irow = it.index();
717 X(irow, j) -= X(jcol, j) * it.value();
723 const MatrixLType& m_mapL;
724 const MatrixUType& m_mapU;
729 template<
typename _MatrixType,
typename Derived,
typename Rhs>
730 struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
731 : solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
733 typedef SparseLU<_MatrixType,Derived> Dec;
734 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
736 template<typename Dest>
void evalTo(Dest& dst)
const
738 dec()._solve(rhs(),dst);
742 template<
typename _MatrixType,
typename Derived,
typename Rhs>
743 struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
744 : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
746 typedef SparseLU<_MatrixType,Derived> Dec;
747 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
749 template<typename Dest>
void evalTo(Dest& dst)
const
751 this->defaultEvalTo(dst);
Index rows() const
Definition: SparseMatrix.h:119
Index cols() const
Definition: SparseMatrix.h:121
const IndicesType & indices() const
Definition: PermutationMatrix.h:358
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, Index > > matrixU() const
Definition: SparseLU.h:144
Transpose< PermutationBase > inverse() const
Definition: PermutationMatrix.h:201
Definition: Constants.h:378
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:207
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:134
Scalar signDeterminant()
Definition: SparseLU.h:311
const internal::solve_retval< SparseLU, Rhs > solve(const MatrixBase< Rhs > &B) const
Definition: SparseLU.h:178
ColXpr col(Index i)
Definition: DenseBase.h:733
std::string lastErrorMessage() const
Definition: SparseLU.h:216
int coletree(const MatrixType &mat, IndexVector &parent, IndexVector &firstRowElt, typename MatrixType::Index *perm=0)
Definition: SparseColEtree.h:61
const PermutationType & colsPermutation() const
Definition: SparseLU.h:161
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:153
Base class of any sparse matrices or sparse expressions.
Definition: SparseMatrixBase.h:26
Derived & derived()
Definition: EigenBase.h:34
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:454
Scalar logAbsDeterminant() const
Definition: SparseLU.h:288
Derived & setConstant(Index size, const Scalar &value)
Definition: CwiseNullaryOp.h:348
void compute(const MatrixType &matrix)
Definition: SparseLU.h:112
Definition: Eigen_Colamd.h:54
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
void isSymmetric(bool sym)
Definition: SparseLU.h:123
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:166
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:371
Expression of a fixed-size or dynamic-size sub-vector.
Definition: ForwardDeclarations.h:83
Definition: Constants.h:376
const unsigned int RowMajorBit
Definition: Constants.h:53
void resize(Index nbRows, Index nbCols)
Definition: PlainObjectBase.h:235
Scalar absDeterminant()
Definition: SparseLU.h:258
Index rows() const
Definition: SparseMatrixBase.h:150
ComputationInfo
Definition: Constants.h:374
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
const internal::sparse_solve_retval< SparseLU, Rhs > solve(const SparseMatrixBase< Rhs > &B) const
Definition: SparseLU.h:191
void treePostorder(Index n, IndexVector &parent, IndexVector &post)
Post order a tree.
Definition: SparseColEtree.h:178
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:515