12#ifndef EIGEN_SPARSE_LU_H
13#define EIGEN_SPARSE_LU_H
17template <
typename _MatrixType,
typename _OrderingType = COLAMDOrdering<
typename _MatrixType::StorageIndex> >
class SparseLU;
18template <
typename MappedSparseMatrixType>
struct SparseLUMatrixLReturnType;
19template <
typename MatrixLType,
typename MatrixUType>
struct SparseLUMatrixUReturnType;
21template <
bool Conjugate,
class SparseLUType>
22class SparseLUTransposeView :
public SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> >
25 typedef SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> > APIBase;
26 using APIBase::m_isInitialized;
28 typedef typename SparseLUType::Scalar Scalar;
29 typedef typename SparseLUType::StorageIndex StorageIndex;
30 typedef typename SparseLUType::MatrixType MatrixType;
31 typedef typename SparseLUType::OrderingType OrderingType;
34 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
35 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
38 SparseLUTransposeView() : m_sparseLU(NULL) {}
39 SparseLUTransposeView(
const SparseLUTransposeView& view) {
40 this->m_sparseLU = view.m_sparseLU;
42 void setIsInitialized(
const bool isInitialized) {this->m_isInitialized = isInitialized;}
43 void setSparseLU(SparseLUType* sparseLU) {m_sparseLU = sparseLU;}
44 using APIBase::_solve_impl;
45 template<
typename Rhs,
typename Dest>
46 bool _solve_impl(
const MatrixBase<Rhs> &B, MatrixBase<Dest> &X_base)
const
48 Dest& X(X_base.derived());
49 eigen_assert(m_sparseLU->info() ==
Success &&
"The matrix should be factorized first");
51 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
55 for(
Index j = 0; j < B.cols(); ++j){
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);
65 for (
Index j = 0; j < B.cols(); ++j)
66 X.col(j) = m_sparseLU->rowsPermutation().transpose() * X.col(j);
69 inline Index rows()
const {
return m_sparseLU->rows(); }
70 inline Index cols()
const {
return m_sparseLU->cols(); }
73 SparseLUType *m_sparseLU;
74 SparseLUTransposeView& operator=(
const SparseLUTransposeView&);
130template <
typename _MatrixType,
typename _OrderingType>
131class SparseLU :
public SparseSolverBase<SparseLU<_MatrixType,_OrderingType> >,
public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::StorageIndex>
135 using APIBase::m_isInitialized;
137 using APIBase::_solve_impl;
139 typedef _MatrixType MatrixType;
140 typedef _OrderingType OrderingType;
141 typedef typename MatrixType::Scalar Scalar;
142 typedef typename MatrixType::RealScalar RealScalar;
143 typedef typename MatrixType::StorageIndex StorageIndex;
145 typedef internal::MappedSuperNodalMatrix<Scalar, StorageIndex> SCMatrix;
149 typedef internal::SparseLUImpl<Scalar, StorageIndex> Base;
152 ColsAtCompileTime = MatrixType::ColsAtCompileTime,
153 MaxColsAtCompileTime = MatrixType::MaxColsAtCompileTime
158 SparseLU():m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(
false),m_diagpivotthresh(1.0),m_detPermR(1)
162 explicit SparseLU(
const MatrixType& matrix)
163 : m_lastError(
""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(
false),m_diagpivotthresh(1.0),m_detPermR(1)
175 void factorize (
const MatrixType& matrix);
176 void simplicialfactorize(
const MatrixType& matrix);
200 const SparseLUTransposeView<false,SparseLU<_MatrixType,_OrderingType> >
transpose()
202 SparseLUTransposeView<false, SparseLU<_MatrixType,_OrderingType> > transposeView;
203 transposeView.setSparseLU(
this);
204 transposeView.setIsInitialized(this->m_isInitialized);
205 return transposeView;
221 const SparseLUTransposeView<true, SparseLU<_MatrixType,_OrderingType> >
adjoint()
223 SparseLUTransposeView<true, SparseLU<_MatrixType,_OrderingType> > adjointView;
224 adjointView.setSparseLU(
this);
225 adjointView.setIsInitialized(this->m_isInitialized);
229 inline Index rows()
const {
return m_mat.
rows(); }
230 inline Index cols()
const {
return m_mat.
cols(); }
234 m_symmetricmode = sym;
243 SparseLUMatrixLReturnType<SCMatrix>
matrixL()
const
245 return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
253 SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >
matrixU()
const
255 return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,StorageIndex> >(m_Lstore, m_Ustore);
277 m_diagpivotthresh = thresh;
280#ifdef EIGEN_PARSED_BY_DOXYGEN
287 template<
typename Rhs>
301 eigen_assert(m_isInitialized &&
"Decomposition is not initialized.");
313 template<
typename Rhs,
typename Dest>
317 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first");
319 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
330 this->
matrixL().solveInPlace(X);
331 this->
matrixU().solveInPlace(X);
353 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
355 Scalar det = Scalar(1.);
358 for (
Index j = 0; j < this->cols(); ++j)
360 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
364 det *= abs(it.value());
385 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
386 Scalar det = Scalar(0.);
387 for (
Index j = 0; j < this->cols(); ++j)
389 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
391 if(it.row() < j)
continue;
394 det += log(abs(it.value()));
408 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
413 for (
Index j = 0; j < this->cols(); ++j)
415 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
421 else if(it.value()==0)
427 return det * m_detPermR * m_detPermC;
436 eigen_assert(m_factorizationIsOk &&
"The matrix should be factorized first.");
438 Scalar det = Scalar(1.);
441 for (
Index j = 0; j < this->cols(); ++j)
443 for (
typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
452 return (m_detPermR * m_detPermC) > 0 ? det : -det;
455 Index nnzL()
const {
return m_nnzL; };
456 Index nnzU()
const {
return m_nnzU; };
460 void initperfvalues()
462 m_perfv.panel_size = 16;
464 m_perfv.maxsuper = 128;
467 m_perfv.fillfactor = 20;
472 bool m_factorizationIsOk;
474 std::string m_lastError;
477 MappedSparseMatrix<Scalar,ColMajor,StorageIndex> m_Ustore;
478 PermutationType m_perm_c;
479 PermutationType m_perm_r ;
482 typename Base::GlobalLU_t m_glu;
485 bool m_symmetricmode;
487 internal::perfvalues m_perfv;
488 RealScalar m_diagpivotthresh;
489 Index m_nnzL, m_nnzU;
490 Index m_detPermR, m_detPermC;
493 SparseLU (
const SparseLU& );
509template <
typename MatrixType,
typename OrderingType>
527 ei_declare_aligned_stack_constructed_variable(StorageIndex,outerIndexPtr,mat.cols()+1,mat.isCompressed()?
const_cast<StorageIndex*
>(mat.outerIndexPtr()):0);
530 if(!mat.isCompressed())
531 IndexVector::Map(outerIndexPtr, mat.cols()+1) = IndexVector::Map(m_mat.outerIndexPtr(),mat.cols()+1);
534 for (
Index i = 0; i < mat.cols(); i++)
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];
543 internal::coletree(m_mat, m_etree,firstRowElt);
546 if (!m_symmetricmode) {
549 internal::treePostorder(StorageIndex(m_mat.cols()), m_etree, post);
553 Index m = m_mat.cols();
555 for (
Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
560 for (
Index i = 0; i < m; i++)
561 post_perm.
indices()(i) = post(i);
564 if(m_perm_c.size()) {
565 m_perm_c = post_perm * m_perm_c;
570 m_analysisIsOk =
true;
594template <
typename MatrixType,
typename OrderingType>
597 using internal::emptyIdxLU;
598 eigen_assert(m_analysisIsOk &&
"analyzePattern() should be called first");
599 eigen_assert((matrix.rows() == matrix.cols()) &&
"Only for squared matrices");
601 m_isInitialized =
true;
610 const StorageIndex * outerIndexPtr;
611 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
614 StorageIndex* outerIndexPtr_t =
new StorageIndex[matrix.cols()+1];
615 for(
Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
616 outerIndexPtr = outerIndexPtr_t;
618 for (
Index i = 0; i < matrix.cols(); i++)
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());
628 for(StorageIndex i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
631 Index m = m_mat.rows();
632 Index n = m_mat.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;
661 tempv.
setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, m) );
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++)
697 if (relax_end(k) != emptyIdxLU)
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++)
775 repfnz_k(irep) = emptyIdxLU;
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;
798template<
typename MappedSupernodalType>
799struct SparseLUMatrixLReturnType : internal::no_assignment_operator
801 typedef typename MappedSupernodalType::Scalar Scalar;
802 explicit SparseLUMatrixLReturnType(
const MappedSupernodalType& mapL) : m_mapL(mapL)
804 Index rows()
const {
return m_mapL.rows(); }
805 Index cols()
const {
return m_mapL.cols(); }
806 template<
typename Dest>
807 void solveInPlace( MatrixBase<Dest> &X)
const
809 m_mapL.solveInPlace(X);
811 template<
bool Conjugate,
typename Dest>
812 void solveTransposedInPlace( MatrixBase<Dest> &X)
const
814 m_mapL.template solveTransposedInPlace<Conjugate>(X);
817 const MappedSupernodalType& m_mapL;
820template<
typename MatrixLType,
typename MatrixUType>
821struct SparseLUMatrixUReturnType : internal::no_assignment_operator
823 typedef typename MatrixLType::Scalar Scalar;
824 SparseLUMatrixUReturnType(
const MatrixLType& mapL,
const MatrixUType& mapU)
825 : m_mapL(mapL),m_mapU(mapU)
827 Index rows()
const {
return m_mapL.rows(); }
828 Index cols()
const {
return m_mapL.cols(); }
830 template<
typename Dest>
void solveInPlace(MatrixBase<Dest> &X)
const
832 Index nrhs = X.cols();
835 for (
Index k = m_mapL.nsuper(); k >= 0; k--)
837 Index fsupc = m_mapL.supToCol()[k];
838 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc];
839 Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
840 Index luptr = m_mapL.colIndexPtr()[fsupc];
844 for (
Index j = 0; j < nrhs; j++)
846 X(fsupc, j) /= m_mapL.valuePtr()[luptr];
852 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
853 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X.coeffRef(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
854 U = A.template triangularView<Upper>().solve(U);
857 for (
Index j = 0; j < nrhs; ++j)
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();
872 template<
bool Conjugate,
typename Dest>
void solveTransposedInPlace(MatrixBase<Dest> &X)
const
875 Index nrhs = X.cols();
878 for (
Index k = 0; k <= m_mapL.nsuper(); k++)
880 Index fsupc = m_mapL.supToCol()[k];
881 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc];
882 Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
883 Index luptr = m_mapL.colIndexPtr()[fsupc];
885 for (
Index j = 0; j < nrhs; ++j)
887 for (
Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
889 typename MatrixUType::InnerIterator it(m_mapU, jcol);
892 Index irow = it.index();
893 X(jcol, j) -= X(irow, j) * (Conjugate?
conj(it.value()): it.value());
899 for (
Index j = 0; j < nrhs; j++)
901 X(fsupc, j) /= (Conjugate?
conj(m_mapL.valuePtr()[luptr]) : m_mapL.valuePtr()[luptr]);
906 Map<const Matrix<Scalar,Dynamic,Dynamic, ColMajor>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
907 Map< Matrix<Scalar,Dynamic,Dest::ColsAtCompileTime, ColMajor>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
909 U = A.adjoint().template triangularView<Lower>().solve(U);
911 U = A.transpose().template triangularView<Lower>().solve(U);
917 const MatrixLType& m_mapL;
918 const MatrixUType& m_mapU;
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Derived & derived()
Definition: EigenBase.h:46
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Sparse matrix.
Definition: MappedSparseMatrix.h:34
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:50
The matrix class, also used for vectors and row-vectors.
Definition: Matrix.h:180
InverseReturnType inverse() const
Definition: PermutationMatrix.h:185
Permutation matrix.
Definition: PermutationMatrix.h:298
const IndicesType & indices() const
Definition: PermutationMatrix.h:360
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
Pseudo expression representing a solving operation.
Definition: Solve.h:63
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:132
Scalar determinant()
Definition: SparseLU.h:434
Scalar absDeterminant()
Definition: SparseLU.h:350
const PermutationType & colsPermutation() const
Definition: SparseLU.h:270
const SparseLUTransposeView< false, SparseLU< _MatrixType, _OrderingType > > transpose()
Definition: SparseLU.h:200
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:595
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:243
std::string lastErrorMessage() const
Definition: SparseLU.h:308
Scalar signDeterminant()
Definition: SparseLU.h:406
const Solve< SparseLU, Rhs > solve(const MatrixBase< Rhs > &B) const
Scalar logAbsDeterminant() const
Definition: SparseLU.h:380
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:275
void compute(const MatrixType &matrix)
Definition: SparseLU.h:182
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:510
const SparseLUTransposeView< true, SparseLU< _MatrixType, _OrderingType > > adjoint()
Definition: SparseLU.h:221
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:299
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:262
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, StorageIndex > > matrixU() const
Definition: SparseLU.h:253
void isSymmetric(bool sym)
Definition: SparseLU.h:232
A versatible sparse matrix representation.
Definition: SparseMatrix.h:98
Index rows() const
Definition: SparseMatrix.h:138
Index cols() const
Definition: SparseMatrix.h:140
A base class for sparse solvers.
Definition: SparseSolverBase.h:68
Expression of a fixed-size or dynamic-size sub-vector.
Definition: VectorBlock.h:60
ComputationInfo
Definition: Constants.h:440
@ NumericalIssue
Definition: Constants.h:444
@ Success
Definition: Constants.h:442
const unsigned int RowMajorBit
Definition: Constants.h:66
Namespace containing all symbols from the Eigen library.
Definition: Core:141
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74