00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef EIGEN_SPARSE_LU_H
00013 #define EIGEN_SPARSE_LU_H
00014
00015 namespace Eigen {
00016
00017 template <typename _MatrixType, typename _OrderingType = COLAMDOrdering<typename _MatrixType::Index> > class SparseLU;
00018 template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType;
00019 template <typename MatrixLType, typename MatrixUType> struct SparseLUMatrixUReturnType;
00020
00072 template <typename _MatrixType, typename _OrderingType>
00073 class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index>
00074 {
00075 public:
00076 typedef _MatrixType MatrixType;
00077 typedef _OrderingType OrderingType;
00078 typedef typename MatrixType::Scalar Scalar;
00079 typedef typename MatrixType::RealScalar RealScalar;
00080 typedef typename MatrixType::Index Index;
00081 typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix;
00082 typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix;
00083 typedef Matrix<Scalar,Dynamic,1> ScalarVector;
00084 typedef Matrix<Index,Dynamic,1> IndexVector;
00085 typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
00086 typedef internal::SparseLUImpl<Scalar, Index> Base;
00087
00088 public:
00089 SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0),m_detPermR(1)
00090 {
00091 initperfvalues();
00092 }
00093 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)
00094 {
00095 initperfvalues();
00096 compute(matrix);
00097 }
00098
00099 ~SparseLU()
00100 {
00101
00102 }
00103
00104 void analyzePattern (const MatrixType& matrix);
00105 void factorize (const MatrixType& matrix);
00106 void simplicialfactorize(const MatrixType& matrix);
00107
00112 void compute (const MatrixType& matrix)
00113 {
00114
00115 analyzePattern(matrix);
00116
00117 factorize(matrix);
00118 }
00119
00120 inline Index rows() const { return m_mat.rows(); }
00121 inline Index cols() const { return m_mat.cols(); }
00123 void isSymmetric(bool sym)
00124 {
00125 m_symmetricmode = sym;
00126 }
00127
00134 SparseLUMatrixLReturnType<SCMatrix> matrixL() const
00135 {
00136 return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
00137 }
00144 SparseLUMatrixUReturnType<SCMatrix,MappedSparseMatrix<Scalar,ColMajor,Index> > matrixU() const
00145 {
00146 return SparseLUMatrixUReturnType<SCMatrix, MappedSparseMatrix<Scalar,ColMajor,Index> >(m_Lstore, m_Ustore);
00147 }
00148
00153 inline const PermutationType& rowsPermutation() const
00154 {
00155 return m_perm_r;
00156 }
00161 inline const PermutationType& colsPermutation() const
00162 {
00163 return m_perm_c;
00164 }
00166 void setPivotThreshold(const RealScalar& thresh)
00167 {
00168 m_diagpivotthresh = thresh;
00169 }
00170
00177 template<typename Rhs>
00178 inline const internal::solve_retval<SparseLU, Rhs> solve(const MatrixBase<Rhs>& B) const
00179 {
00180 eigen_assert(m_factorizationIsOk && "SparseLU is not initialized.");
00181 eigen_assert(rows()==B.rows()
00182 && "SparseLU::solve(): invalid number of rows of the right hand side matrix B");
00183 return internal::solve_retval<SparseLU, Rhs>(*this, B.derived());
00184 }
00185
00190 template<typename Rhs>
00191 inline const internal::sparse_solve_retval<SparseLU, Rhs> solve(const SparseMatrixBase<Rhs>& B) const
00192 {
00193 eigen_assert(m_factorizationIsOk && "SparseLU is not initialized.");
00194 eigen_assert(rows()==B.rows()
00195 && "SparseLU::solve(): invalid number of rows of the right hand side matrix B");
00196 return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived());
00197 }
00198
00207 ComputationInfo info() const
00208 {
00209 eigen_assert(m_isInitialized && "Decomposition is not initialized.");
00210 return m_info;
00211 }
00212
00216 std::string lastErrorMessage() const
00217 {
00218 return m_lastError;
00219 }
00220
00221 template<typename Rhs, typename Dest>
00222 bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &_X) const
00223 {
00224 Dest& X(_X.derived());
00225 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first");
00226 EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
00227 THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
00228
00229
00230
00231 X.resize(B.rows(),B.cols());
00232 for(Index j = 0; j < B.cols(); ++j)
00233 X.col(j) = rowsPermutation() * B.col(j);
00234
00235
00236 this->matrixL().solveInPlace(X);
00237 this->matrixU().solveInPlace(X);
00238
00239
00240 for (Index j = 0; j < B.cols(); ++j)
00241 X.col(j) = colsPermutation().inverse() * X.col(j);
00242
00243 return true;
00244 }
00245
00256 Scalar absDeterminant()
00257 {
00258 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
00259
00260 Scalar det = Scalar(1.);
00261
00262
00263 for (Index j = 0; j < this->cols(); ++j)
00264 {
00265 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
00266 {
00267 if(it.row() < j) continue;
00268 if(it.row() == j)
00269 {
00270 det *= (std::abs)(it.value());
00271 break;
00272 }
00273 }
00274 }
00275 return det;
00276 }
00277
00286 Scalar logAbsDeterminant() const
00287 {
00288 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
00289 Scalar det = Scalar(0.);
00290 for (Index j = 0; j < this->cols(); ++j)
00291 {
00292 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
00293 {
00294 if(it.row() < j) continue;
00295 if(it.row() == j)
00296 {
00297 det += (std::log)((std::abs)(it.value()));
00298 break;
00299 }
00300 }
00301 }
00302 return det;
00303 }
00304
00309 Scalar signDeterminant()
00310 {
00311 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
00312 return Scalar(m_detPermR);
00313 }
00314
00315 protected:
00316
00317 void initperfvalues()
00318 {
00319 m_perfv.panel_size = 1;
00320 m_perfv.relax = 1;
00321 m_perfv.maxsuper = 128;
00322 m_perfv.rowblk = 16;
00323 m_perfv.colblk = 8;
00324 m_perfv.fillfactor = 20;
00325 }
00326
00327
00328 mutable ComputationInfo m_info;
00329 bool m_isInitialized;
00330 bool m_factorizationIsOk;
00331 bool m_analysisIsOk;
00332 std::string m_lastError;
00333 NCMatrix m_mat;
00334 SCMatrix m_Lstore;
00335 MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore;
00336 PermutationType m_perm_c;
00337 PermutationType m_perm_r ;
00338 IndexVector m_etree;
00339
00340 typename Base::GlobalLU_t m_glu;
00341
00342
00343 bool m_symmetricmode;
00344
00345 internal::perfvalues<Index> m_perfv;
00346 RealScalar m_diagpivotthresh;
00347 Index m_nnzL, m_nnzU;
00348 Index m_detPermR;
00349 private:
00350
00351 SparseLU (const SparseLU& );
00352
00353 };
00354
00355
00356
00357
00368 template <typename MatrixType, typename OrderingType>
00369 void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
00370 {
00371
00372
00373
00374 OrderingType ord;
00375 ord(mat,m_perm_c);
00376
00377
00378
00379 m_mat = mat;
00380 if (m_perm_c.size()) {
00381 m_mat.uncompress();
00382
00383 const Index * outerIndexPtr;
00384 if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr();
00385 else
00386 {
00387 Index *outerIndexPtr_t = new Index[mat.cols()+1];
00388 for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
00389 outerIndexPtr = outerIndexPtr_t;
00390 }
00391 for (Index i = 0; i < mat.cols(); i++)
00392 {
00393 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
00394 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
00395 }
00396 if(!mat.isCompressed()) delete[] outerIndexPtr;
00397 }
00398
00399 IndexVector firstRowElt;
00400 internal::coletree(m_mat, m_etree,firstRowElt);
00401
00402
00403 if (!m_symmetricmode) {
00404 IndexVector post, iwork;
00405
00406 internal::treePostorder(m_mat.cols(), m_etree, post);
00407
00408
00409
00410 Index m = m_mat.cols();
00411 iwork.resize(m+1);
00412 for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i));
00413 m_etree = iwork;
00414
00415
00416 PermutationType post_perm(m);
00417 for (Index i = 0; i < m; i++)
00418 post_perm.indices()(i) = post(i);
00419
00420
00421 if(m_perm_c.size()) {
00422 m_perm_c = post_perm * m_perm_c;
00423 }
00424
00425 }
00426
00427 m_analysisIsOk = true;
00428 }
00429
00430
00431
00432
00451 template <typename MatrixType, typename OrderingType>
00452 void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
00453 {
00454 using internal::emptyIdxLU;
00455 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first");
00456 eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
00457
00458 typedef typename IndexVector::Scalar Index;
00459
00460
00461
00462
00463 m_mat = matrix;
00464 if (m_perm_c.size())
00465 {
00466 m_mat.uncompress();
00467
00468 const Index * outerIndexPtr;
00469 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr();
00470 else
00471 {
00472 Index* outerIndexPtr_t = new Index[matrix.cols()+1];
00473 for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i];
00474 outerIndexPtr = outerIndexPtr_t;
00475 }
00476 for (Index i = 0; i < matrix.cols(); i++)
00477 {
00478 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i];
00479 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i];
00480 }
00481 if(!matrix.isCompressed()) delete[] outerIndexPtr;
00482 }
00483 else
00484 {
00485 m_perm_c.resize(matrix.cols());
00486 for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
00487 }
00488
00489 Index m = m_mat.rows();
00490 Index n = m_mat.cols();
00491 Index nnz = m_mat.nonZeros();
00492 Index maxpanel = m_perfv.panel_size * m;
00493
00494 Index lwork = 0;
00495 Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu);
00496 if (info)
00497 {
00498 m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
00499 m_factorizationIsOk = false;
00500 return ;
00501 }
00502
00503
00504 IndexVector segrep(m); segrep.setZero();
00505 IndexVector parent(m); parent.setZero();
00506 IndexVector xplore(m); xplore.setZero();
00507 IndexVector repfnz(maxpanel);
00508 IndexVector panel_lsub(maxpanel);
00509 IndexVector xprune(n); xprune.setZero();
00510 IndexVector marker(m*internal::LUNoMarker); marker.setZero();
00511
00512 repfnz.setConstant(-1);
00513 panel_lsub.setConstant(-1);
00514
00515
00516 ScalarVector dense;
00517 dense.setZero(maxpanel);
00518 ScalarVector tempv;
00519 tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, m) );
00520
00521
00522 PermutationType iperm_c(m_perm_c.inverse());
00523
00524
00525 IndexVector relax_end(n);
00526 if ( m_symmetricmode == true )
00527 Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
00528 else
00529 Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
00530
00531
00532 m_perm_r.resize(m);
00533 m_perm_r.indices().setConstant(-1);
00534 marker.setConstant(-1);
00535 m_detPermR = 1;
00536
00537 m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
00538 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
00539
00540
00541
00542
00543 Index jcol;
00544 IndexVector panel_histo(n);
00545 Index pivrow;
00546 Index nseg1;
00547 Index nseg;
00548 Index irep;
00549 Index i, k, jj;
00550 for (jcol = 0; jcol < n; )
00551 {
00552
00553 Index panel_size = m_perfv.panel_size;
00554 for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
00555 {
00556 if (relax_end(k) != emptyIdxLU)
00557 {
00558 panel_size = k - jcol;
00559 break;
00560 }
00561 }
00562 if (k == n)
00563 panel_size = n - jcol;
00564
00565
00566 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);
00567
00568
00569 Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu);
00570
00571
00572 for ( jj = jcol; jj< jcol + panel_size; jj++)
00573 {
00574 k = (jj - jcol) * m;
00575
00576 nseg = nseg1;
00577
00578 VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
00579 VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
00580 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);
00581 if ( info )
00582 {
00583 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
00584 m_info = NumericalIssue;
00585 m_factorizationIsOk = false;
00586 return;
00587 }
00588
00589 VectorBlock<ScalarVector> dense_k(dense, k, m);
00590 VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1);
00591 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu);
00592 if ( info )
00593 {
00594 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
00595 m_info = NumericalIssue;
00596 m_factorizationIsOk = false;
00597 return;
00598 }
00599
00600
00601 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu);
00602 if ( info )
00603 {
00604 m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
00605 m_info = NumericalIssue;
00606 m_factorizationIsOk = false;
00607 return;
00608 }
00609
00610
00611 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
00612 if ( info )
00613 {
00614 m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
00615 std::ostringstream returnInfo;
00616 returnInfo << info;
00617 m_lastError += returnInfo.str();
00618 m_info = NumericalIssue;
00619 m_factorizationIsOk = false;
00620 return;
00621 }
00622
00623
00624 if (pivrow != jj) m_detPermR *= -1;
00625
00626
00627 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu);
00628
00629
00630 for (i = 0; i < nseg; i++)
00631 {
00632 irep = segrep(i);
00633 repfnz_k(irep) = emptyIdxLU;
00634 }
00635 }
00636 jcol += panel_size;
00637 }
00638
00639
00640 Base::countnz(n, m_nnzL, m_nnzU, m_glu);
00641
00642 Base::fixupL(n, m_perm_r.indices(), m_glu);
00643
00644
00645 m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup);
00646
00647 new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, Index> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() );
00648
00649 m_info = Success;
00650 m_factorizationIsOk = true;
00651 }
00652
00653 template<typename MappedSupernodalType>
00654 struct SparseLUMatrixLReturnType : internal::no_assignment_operator
00655 {
00656 typedef typename MappedSupernodalType::Index Index;
00657 typedef typename MappedSupernodalType::Scalar Scalar;
00658 SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL)
00659 { }
00660 Index rows() { return m_mapL.rows(); }
00661 Index cols() { return m_mapL.cols(); }
00662 template<typename Dest>
00663 void solveInPlace( MatrixBase<Dest> &X) const
00664 {
00665 m_mapL.solveInPlace(X);
00666 }
00667 const MappedSupernodalType& m_mapL;
00668 };
00669
00670 template<typename MatrixLType, typename MatrixUType>
00671 struct SparseLUMatrixUReturnType : internal::no_assignment_operator
00672 {
00673 typedef typename MatrixLType::Index Index;
00674 typedef typename MatrixLType::Scalar Scalar;
00675 SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU)
00676 : m_mapL(mapL),m_mapU(mapU)
00677 { }
00678 Index rows() { return m_mapL.rows(); }
00679 Index cols() { return m_mapL.cols(); }
00680
00681 template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const
00682 {
00683 Index nrhs = X.cols();
00684 Index n = X.rows();
00685
00686 for (Index k = m_mapL.nsuper(); k >= 0; k--)
00687 {
00688 Index fsupc = m_mapL.supToCol()[k];
00689 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc];
00690 Index nsupc = m_mapL.supToCol()[k+1] - fsupc;
00691 Index luptr = m_mapL.colIndexPtr()[fsupc];
00692
00693 if (nsupc == 1)
00694 {
00695 for (Index j = 0; j < nrhs; j++)
00696 {
00697 X(fsupc, j) /= m_mapL.valuePtr()[luptr];
00698 }
00699 }
00700 else
00701 {
00702 Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) );
00703 Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) );
00704 U = A.template triangularView<Upper>().solve(U);
00705 }
00706
00707 for (Index j = 0; j < nrhs; ++j)
00708 {
00709 for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++)
00710 {
00711 typename MatrixUType::InnerIterator it(m_mapU, jcol);
00712 for ( ; it; ++it)
00713 {
00714 Index irow = it.index();
00715 X(irow, j) -= X(jcol, j) * it.value();
00716 }
00717 }
00718 }
00719 }
00720 }
00721 const MatrixLType& m_mapL;
00722 const MatrixUType& m_mapU;
00723 };
00724
00725 namespace internal {
00726
00727 template<typename _MatrixType, typename Derived, typename Rhs>
00728 struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
00729 : solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
00730 {
00731 typedef SparseLU<_MatrixType,Derived> Dec;
00732 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs)
00733
00734 template<typename Dest> void evalTo(Dest& dst) const
00735 {
00736 dec()._solve(rhs(),dst);
00737 }
00738 };
00739
00740 template<typename _MatrixType, typename Derived, typename Rhs>
00741 struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
00742 : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
00743 {
00744 typedef SparseLU<_MatrixType,Derived> Dec;
00745 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
00746
00747 template<typename Dest> void evalTo(Dest& dst) const
00748 {
00749 this->defaultEvalTo(dst);
00750 }
00751 };
00752 }
00753
00754 }
00755
00756 #endif