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;
21 template <
bool Conjugate,
class SparseLUType>
22 class SparseLUTransposeView :
public SparseSolverBase<SparseLUTransposeView<Conjugate,SparseLUType> >
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>
49 eigen_assert(m_sparseLU->info() ==
Success &&
"The matrix should be factorized first");
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);
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&);
130 template <
typename _MatrixType,
typename _OrderingType>
131 class 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)
174 void analyzePattern (
const MatrixType& matrix);
175 void factorize (
const MatrixType& matrix);
176 void simplicialfactorize(
const MatrixType& matrix);
185 analyzePattern(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>
289 #endif // EIGEN_PARSED_BY_DOXYGEN 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);
327 X.col(j) = rowsPermutation() * B.const_cast_derived().col(j);
330 this->matrixL().solveInPlace(X);
331 this->matrixU().solveInPlace(X);
335 X.col(j) = colsPermutation().inverse() * X.col(j);
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;
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;
509 template <
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;
594 template <
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;
798 template<
typename MappedSupernodalType>
799 struct 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>
809 m_mapL.solveInPlace(X);
811 template<
bool Conjugate,
typename Dest>
814 m_mapL.template solveTransposedInPlace<Conjugate>(X);
817 const MappedSupernodalType& m_mapL;
820 template<
typename MatrixLType,
typename MatrixUType>
821 struct 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(); }
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];
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 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]);
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;
const PermutationType & rowsPermutation() const
Definition: SparseLU.h:262
A matrix or vector expression mapping an existing array of data.
Definition: Map.h:94
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_conjugate_op< typename Derived::Scalar >, const Derived > conj(const Eigen::ArrayBase< Derived > &x)
A base class for sparse solvers.
Definition: SparseSolverBase.h:67
Namespace containing all symbols from the Eigen library.
Definition: Core:141
void factorize(const MatrixType &matrix)
Definition: SparseLU.h:595
Derived & derived()
Definition: EigenBase.h:46
const unsigned int RowMajorBit
Definition: Constants.h:66
void resize(Index rows, Index cols)
Definition: PlainObjectBase.h:271
Sparse supernodal LU factorization for general matrices.
Definition: SparseLU.h:17
void setPivotThreshold(const RealScalar &thresh)
Definition: SparseLU.h:275
Scalar absDeterminant()
Definition: SparseLU.h:350
Derived & setZero(Index size)
Definition: CwiseNullaryOp.h:562
SparseLUMatrixLReturnType< SCMatrix > matrixL() const
Definition: SparseLU.h:243
Expression of a fixed-size or dynamic-size sub-vector.
Definition: VectorBlock.h:56
const SparseLUTransposeView< true, SparseLU< _MatrixType, _OrderingType > > adjoint()
Definition: SparseLU.h:221
const IndicesType & indices() const
Definition: PermutationMatrix.h:360
void compute(const MatrixType &matrix)
Definition: SparseLU.h:182
void isSymmetric(bool sym)
Definition: SparseLU.h:232
Definition: Constants.h:444
EIGEN_DEFAULT_DENSE_INDEX_TYPE Index
The Index type as used for the API.
Definition: Meta.h:74
EIGEN_CONSTEXPR Index cols() const EIGEN_NOEXCEPT
Definition: EigenBase.h:63
Scalar determinant()
Definition: SparseLU.h:434
Derived & setConstant(Index size, const Scalar &val)
Definition: CwiseNullaryOp.h:361
Definition: Constants.h:442
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_abs_op< typename Derived::Scalar >, const Derived > abs(const Eigen::ArrayBase< Derived > &x)
SparseLUMatrixUReturnType< SCMatrix, MappedSparseMatrix< Scalar, ColMajor, StorageIndex > > matrixU() const
Definition: SparseLU.h:253
const Eigen::CwiseUnaryOp< Eigen::internal::scalar_log_op< typename Derived::Scalar >, const Derived > log(const Eigen::ArrayBase< Derived > &x)
std::string lastErrorMessage() const
Definition: SparseLU.h:308
Scalar logAbsDeterminant() const
Definition: SparseLU.h:380
const SparseLUTransposeView< false, SparseLU< _MatrixType, _OrderingType > > transpose()
Definition: SparseLU.h:200
const PermutationType & colsPermutation() const
Definition: SparseLU.h:270
CoeffReturnType value() const
Definition: DenseBase.h:526
EIGEN_CONSTEXPR Index rows() const EIGEN_NOEXCEPT
Definition: EigenBase.h:60
Pseudo expression representing a solving operation.
Definition: Solve.h:62
void analyzePattern(const MatrixType &matrix)
Definition: SparseLU.h:510
Convenience specialization of Stride to specify only an outer stride See class Map for some examples...
Definition: Stride.h:106
ComputationInfo
Definition: Constants.h:440
Scalar & coeffRef(Index row, Index col)
Definition: DenseCoeffsBase.h:342
Scalar signDeterminant()
Definition: SparseLU.h:406
Base class for all dense matrices, vectors, and expressions.
Definition: MatrixBase.h:48
ComputationInfo info() const
Reports whether previous computation was successful.
Definition: SparseLU.h:299