12 #ifndef EIGEN_SPARSE_LU_H
13 #define EIGEN_SPARSE_LU_H
17 template <
typename _MatrixType,
typename _OrderingType = COLAMDOrdering<
typename _MatrixType::StorageIndex> >
class SparseLU;
21 template <
bool Conjugate,
class SparseLUType>
45 template<
typename Rhs,
typename Dest>
48 Dest&
X(X_base.derived());
51 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
56 X.col(
j) =
m_sparseLU->colsPermutation() *
B.const_cast_derived().col(
j);
59 m_sparseLU->matrixU().template solveTransposedInPlace<Conjugate>(
X);
62 m_sparseLU->matrixL().template solveTransposedInPlace<Conjugate>(
X);
130 template <
typename _MatrixType,
typename _OrderingType>
205 return transposeView;
280 #ifdef EIGEN_PARSED_BY_DOXYGEN
287 template<
typename Rhs>
289 #endif // EIGEN_PARSED_BY_DOXYGEN
301 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
313 template<
typename Rhs,
typename Dest>
316 Dest&
X(X_base.derived());
319 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
323 X.resize(
B.rows(),
B.cols());
330 this->
matrixL().solveInPlace(X);
331 this->
matrixU().solveInPlace(X);
360 for (
typename SCMatrix::InnerIterator it(
m_Lstore,
j); it; ++it)
364 det *=
abs(it.value());
389 for (
typename SCMatrix::InnerIterator it(
m_Lstore,
j); it; ++it)
391 if(it.row() <
j)
continue;
394 det +=
log(
abs(it.value()));
415 for (
typename SCMatrix::InnerIterator it(
m_Lstore,
j); it; ++it)
421 else if(it.value()==0)
443 for (
typename SCMatrix::InnerIterator it(
m_Lstore,
j); it; ++it)
509 template <
typename MatrixType,
typename OrderingType>
530 if(!
mat.isCompressed())
531 IndexVector::Map(outerIndexPtr,
mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),
mat.cols()+1);
536 m_mat.outerIndexPtr()[m_perm_c.indices()(
i)] = outerIndexPtr[
i];
537 m_mat.innerNonZeroPtr()[m_perm_c.indices()(
i)] = outerIndexPtr[
i+1] - outerIndexPtr[
i];
546 if (!m_symmetricmode) {
555 for (
Index i = 0;
i <
m; ++
i) iwork(post(
i)) = post(m_etree(
i));
564 if(m_perm_c.size()) {
565 m_perm_c = post_perm * m_perm_c;
570 m_analysisIsOk =
true;
594 template <
typename MatrixType,
typename OrderingType>
598 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
601 m_isInitialized =
true;
611 if (
matrix.isCompressed()) outerIndexPtr =
matrix.outerIndexPtr();
615 for(
Index i = 0;
i <=
matrix.cols();
i++) outerIndexPtr_t[
i] = m_mat.outerIndexPtr()[
i];
616 outerIndexPtr = outerIndexPtr_t;
620 m_mat.outerIndexPtr()[m_perm_c.indices()(
i)] = outerIndexPtr[
i];
621 m_mat.innerNonZeroPtr()[m_perm_c.indices()(
i)] = outerIndexPtr[
i+1] - outerIndexPtr[
i];
623 if(!
matrix.isCompressed())
delete[] outerIndexPtr;
627 m_perm_c.resize(
matrix.cols());
633 Index nnz = m_mat.nonZeros();
634 Index maxpanel = m_perfv.panel_size *
m;
637 Index info = Base::memInit(
m,
n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
640 m_lastError =
"UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
641 m_factorizationIsOk =
false;
668 if ( m_symmetricmode ==
true )
669 Base::heap_relax_snode(
n, m_etree, m_perfv.relax, marker, relax_end);
671 Base::relax_snode(
n, m_etree, m_perfv.relax, marker, relax_end);
675 m_perm_r.indices().setConstant(-1);
679 m_glu.supno(0) =
emptyIdxLU; m_glu.xsup.setConstant(0);
680 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) =
Index(0);
691 for (jcol = 0; jcol <
n; )
694 Index panel_size = m_perfv.panel_size;
695 for (k = jcol + 1; k < (
std::min)(jcol+panel_size,
n); k++)
699 panel_size = k - jcol;
704 panel_size =
n - jcol;
707 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);
710 Base::panel_bmod(
m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
713 for ( jj = jcol; jj< jcol + panel_size; jj++)
721 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);
724 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
726 m_factorizationIsOk =
false;
732 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
735 m_lastError =
"UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
737 m_factorizationIsOk =
false;
742 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
745 m_lastError =
"UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
747 m_factorizationIsOk =
false;
752 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.
indices(), pivrow, m_glu);
755 m_lastError =
"THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
756 std::ostringstream returnInfo;
758 m_lastError += returnInfo.str();
760 m_factorizationIsOk =
false;
766 if (pivrow != jj) m_detPermR = -m_detPermR;
769 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
772 for (
i = 0;
i < nseg;
i++)
781 m_detPermR = m_perm_r.determinant();
782 m_detPermC = m_perm_c.determinant();
785 Base::countnz(
n, m_nnzL, m_nnzU, m_glu);
787 Base::fixupL(
n, m_perm_r.indices(), m_glu);
790 m_Lstore.setInfos(
m,
n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
795 m_factorizationIsOk =
true;
798 template<
typename MappedSupernodalType>
806 template<
typename Dest>
811 template<
bool Conjugate,
typename Dest>
814 m_mapL.template solveTransposedInPlace<Conjugate>(
X);
820 template<
typename MatrixLType,
typename MatrixUType>
846 X(fsupc,
j) /=
m_mapL.valuePtr()[luptr];
854 U =
A.template triangularView<Upper>().solve(
U);
859 for (
Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
861 typename MatrixUType::InnerIterator it(
m_mapU, jcol);
864 Index irow = it.index();
865 X(irow,
j) -=
X(jcol,
j) * it.value();
887 for (
Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
889 typename MatrixUType::InnerIterator it(
m_mapU, jcol);
892 Index irow = it.index();