$treeview $search $mathjax
Eigen
3.2.5
$projectbrief
|
$projectbrief
|
$searchbox |
00001 // This file is part of Eigen, a lightweight C++ template library 00002 // for linear algebra. 00003 // 00004 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr> 00005 // Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr> 00006 // 00007 // This Source Code Form is subject to the terms of the Mozilla 00008 // Public License v. 2.0. If a copy of the MPL was not distributed 00009 // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. 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 // Free all explicit dynamic pointers 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 // Analyze 00115 analyzePattern(matrix); 00116 //Factorize 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_base) const 00223 { 00224 Dest& X(X_base.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 // Permute the right hand side to form X = Pr*B 00230 // on return, X is overwritten by the computed solution 00231 X.resize(B.rows(),B.cols()); 00232 00233 // this ugly const_cast_derived() helps to detect aliasing when applying the permutations 00234 for(Index j = 0; j < B.cols(); ++j) 00235 X.col(j) = rowsPermutation() * B.const_cast_derived().col(j); 00236 00237 //Forward substitution with L 00238 this->matrixL().solveInPlace(X); 00239 this->matrixU().solveInPlace(X); 00240 00241 // Permute back the solution 00242 for (Index j = 0; j < B.cols(); ++j) 00243 X.col(j) = colsPermutation().inverse() * X.col(j); 00244 00245 return true; 00246 } 00247 00258 Scalar absDeterminant() 00259 { 00260 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 00261 // Initialize with the determinant of the row matrix 00262 Scalar det = Scalar(1.); 00263 // Note that the diagonal blocks of U are stored in supernodes, 00264 // which are available in the L part :) 00265 for (Index j = 0; j < this->cols(); ++j) 00266 { 00267 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 00268 { 00269 if(it.index() == j) 00270 { 00271 using std::abs; 00272 det *= abs(it.value()); 00273 break; 00274 } 00275 } 00276 } 00277 return det; 00278 } 00279 00288 Scalar logAbsDeterminant() const 00289 { 00290 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 00291 Scalar det = Scalar(0.); 00292 for (Index j = 0; j < this->cols(); ++j) 00293 { 00294 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 00295 { 00296 if(it.row() < j) continue; 00297 if(it.row() == j) 00298 { 00299 using std::log; using std::abs; 00300 det += log(abs(it.value())); 00301 break; 00302 } 00303 } 00304 } 00305 return det; 00306 } 00307 00312 Scalar signDeterminant() 00313 { 00314 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 00315 // Initialize with the determinant of the row matrix 00316 Index det = 1; 00317 // Note that the diagonal blocks of U are stored in supernodes, 00318 // which are available in the L part :) 00319 for (Index j = 0; j < this->cols(); ++j) 00320 { 00321 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 00322 { 00323 if(it.index() == j) 00324 { 00325 if(it.value()<0) 00326 det = -det; 00327 else if(it.value()==0) 00328 return 0; 00329 break; 00330 } 00331 } 00332 } 00333 return det * m_detPermR * m_detPermC; 00334 } 00335 00340 Scalar determinant() 00341 { 00342 eigen_assert(m_factorizationIsOk && "The matrix should be factorized first."); 00343 // Initialize with the determinant of the row matrix 00344 Scalar det = Scalar(1.); 00345 // Note that the diagonal blocks of U are stored in supernodes, 00346 // which are available in the L part :) 00347 for (Index j = 0; j < this->cols(); ++j) 00348 { 00349 for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it) 00350 { 00351 if(it.index() == j) 00352 { 00353 det *= it.value(); 00354 break; 00355 } 00356 } 00357 } 00358 return det * Scalar(m_detPermR * m_detPermC); 00359 } 00360 00361 protected: 00362 // Functions 00363 void initperfvalues() 00364 { 00365 m_perfv.panel_size = 16; 00366 m_perfv.relax = 1; 00367 m_perfv.maxsuper = 128; 00368 m_perfv.rowblk = 16; 00369 m_perfv.colblk = 8; 00370 m_perfv.fillfactor = 20; 00371 } 00372 00373 // Variables 00374 mutable ComputationInfo m_info; 00375 bool m_isInitialized; 00376 bool m_factorizationIsOk; 00377 bool m_analysisIsOk; 00378 std::string m_lastError; 00379 NCMatrix m_mat; // The input (permuted ) matrix 00380 SCMatrix m_Lstore; // The lower triangular matrix (supernodal) 00381 MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore; // The upper triangular matrix 00382 PermutationType m_perm_c; // Column permutation 00383 PermutationType m_perm_r ; // Row permutation 00384 IndexVector m_etree; // Column elimination tree 00385 00386 typename Base::GlobalLU_t m_glu; 00387 00388 // SparseLU options 00389 bool m_symmetricmode; 00390 // values for performance 00391 internal::perfvalues<Index> m_perfv; 00392 RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot 00393 Index m_nnzL, m_nnzU; // Nonzeros in L and U factors 00394 Index m_detPermR, m_detPermC; // Determinants of the permutation matrices 00395 private: 00396 // Disable copy constructor 00397 SparseLU (const SparseLU& ); 00398 00399 }; // End class SparseLU 00400 00401 00402 00403 // Functions needed by the anaysis phase 00414 template <typename MatrixType, typename OrderingType> 00415 void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat) 00416 { 00417 00418 //TODO It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat. 00419 00420 OrderingType ord; 00421 ord(mat,m_perm_c); 00422 00423 // Apply the permutation to the column of the input matrix 00424 //First copy the whole input matrix. 00425 m_mat = mat; 00426 if (m_perm_c.size()) { 00427 m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used. 00428 //Then, permute only the column pointers 00429 const Index * outerIndexPtr; 00430 if (mat.isCompressed()) outerIndexPtr = mat.outerIndexPtr(); 00431 else 00432 { 00433 Index *outerIndexPtr_t = new Index[mat.cols()+1]; 00434 for(Index i = 0; i <= mat.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 00435 outerIndexPtr = outerIndexPtr_t; 00436 } 00437 for (Index i = 0; i < mat.cols(); i++) 00438 { 00439 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 00440 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 00441 } 00442 if(!mat.isCompressed()) delete[] outerIndexPtr; 00443 } 00444 // Compute the column elimination tree of the permuted matrix 00445 IndexVector firstRowElt; 00446 internal::coletree(m_mat, m_etree,firstRowElt); 00447 00448 // In symmetric mode, do not do postorder here 00449 if (!m_symmetricmode) { 00450 IndexVector post, iwork; 00451 // Post order etree 00452 internal::treePostorder(m_mat.cols(), m_etree, post); 00453 00454 00455 // Renumber etree in postorder 00456 Index m = m_mat.cols(); 00457 iwork.resize(m+1); 00458 for (Index i = 0; i < m; ++i) iwork(post(i)) = post(m_etree(i)); 00459 m_etree = iwork; 00460 00461 // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree 00462 PermutationType post_perm(m); 00463 for (Index i = 0; i < m; i++) 00464 post_perm.indices()(i) = post(i); 00465 00466 // Combine the two permutations : postorder the permutation for future use 00467 if(m_perm_c.size()) { 00468 m_perm_c = post_perm * m_perm_c; 00469 } 00470 00471 } // end postordering 00472 00473 m_analysisIsOk = true; 00474 } 00475 00476 // Functions needed by the numerical factorization phase 00477 00478 00497 template <typename MatrixType, typename OrderingType> 00498 void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix) 00499 { 00500 using internal::emptyIdxLU; 00501 eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 00502 eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices"); 00503 00504 typedef typename IndexVector::Scalar Index; 00505 00506 00507 // Apply the column permutation computed in analyzepattern() 00508 // m_mat = matrix * m_perm_c.inverse(); 00509 m_mat = matrix; 00510 if (m_perm_c.size()) 00511 { 00512 m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. 00513 //Then, permute only the column pointers 00514 const Index * outerIndexPtr; 00515 if (matrix.isCompressed()) outerIndexPtr = matrix.outerIndexPtr(); 00516 else 00517 { 00518 Index* outerIndexPtr_t = new Index[matrix.cols()+1]; 00519 for(Index i = 0; i <= matrix.cols(); i++) outerIndexPtr_t[i] = m_mat.outerIndexPtr()[i]; 00520 outerIndexPtr = outerIndexPtr_t; 00521 } 00522 for (Index i = 0; i < matrix.cols(); i++) 00523 { 00524 m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i]; 00525 m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = outerIndexPtr[i+1] - outerIndexPtr[i]; 00526 } 00527 if(!matrix.isCompressed()) delete[] outerIndexPtr; 00528 } 00529 else 00530 { //FIXME This should not be needed if the empty permutation is handled transparently 00531 m_perm_c.resize(matrix.cols()); 00532 for(Index i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i; 00533 } 00534 00535 Index m = m_mat.rows(); 00536 Index n = m_mat.cols(); 00537 Index nnz = m_mat.nonZeros(); 00538 Index maxpanel = m_perfv.panel_size * m; 00539 // Allocate working storage common to the factor routines 00540 Index lwork = 0; 00541 Index info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 00542 if (info) 00543 { 00544 m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ; 00545 m_factorizationIsOk = false; 00546 return ; 00547 } 00548 00549 // Set up pointers for integer working arrays 00550 IndexVector segrep(m); segrep.setZero(); 00551 IndexVector parent(m); parent.setZero(); 00552 IndexVector xplore(m); xplore.setZero(); 00553 IndexVector repfnz(maxpanel); 00554 IndexVector panel_lsub(maxpanel); 00555 IndexVector xprune(n); xprune.setZero(); 00556 IndexVector marker(m*internal::LUNoMarker); marker.setZero(); 00557 00558 repfnz.setConstant(-1); 00559 panel_lsub.setConstant(-1); 00560 00561 // Set up pointers for scalar working arrays 00562 ScalarVector dense; 00563 dense.setZero(maxpanel); 00564 ScalarVector tempv; 00565 tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) ); 00566 00567 // Compute the inverse of perm_c 00568 PermutationType iperm_c(m_perm_c.inverse()); 00569 00570 // Identify initial relaxed snodes 00571 IndexVector relax_end(n); 00572 if ( m_symmetricmode == true ) 00573 Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 00574 else 00575 Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end); 00576 00577 00578 m_perm_r.resize(m); 00579 m_perm_r.indices().setConstant(-1); 00580 marker.setConstant(-1); 00581 m_detPermR = 1; // Record the determinant of the row permutation 00582 00583 m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0); 00584 m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0); 00585 00586 // Work on one 'panel' at a time. A panel is one of the following : 00587 // (a) a relaxed supernode at the bottom of the etree, or 00588 // (b) panel_size contiguous columns, <panel_size> defined by the user 00589 Index jcol; 00590 IndexVector panel_histo(n); 00591 Index pivrow; // Pivotal row number in the original row matrix 00592 Index nseg1; // Number of segments in U-column above panel row jcol 00593 Index nseg; // Number of segments in each U-column 00594 Index irep; 00595 Index i, k, jj; 00596 for (jcol = 0; jcol < n; ) 00597 { 00598 // Adjust panel size so that a panel won't overlap with the next relaxed snode. 00599 Index panel_size = m_perfv.panel_size; // upper bound on panel width 00600 for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++) 00601 { 00602 if (relax_end(k) != emptyIdxLU) 00603 { 00604 panel_size = k - jcol; 00605 break; 00606 } 00607 } 00608 if (k == n) 00609 panel_size = n - jcol; 00610 00611 // Symbolic outer factorization on a panel of columns 00612 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); 00613 00614 // Numeric sup-panel updates in topological order 00615 Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 00616 00617 // Sparse LU within the panel, and below the panel diagonal 00618 for ( jj = jcol; jj< jcol + panel_size; jj++) 00619 { 00620 k = (jj - jcol) * m; // Column index for w-wide arrays 00621 00622 nseg = nseg1; // begin after all the panel segments 00623 //Depth-first-search for the current column 00624 VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); 00625 VectorBlock<IndexVector> repfnz_k(repfnz, k, m); 00626 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); 00627 if ( info ) 00628 { 00629 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() "; 00630 m_info = NumericalIssue; 00631 m_factorizationIsOk = false; 00632 return; 00633 } 00634 // Numeric updates to this column 00635 VectorBlock<ScalarVector> dense_k(dense, k, m); 00636 VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1); 00637 info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 00638 if ( info ) 00639 { 00640 m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() "; 00641 m_info = NumericalIssue; 00642 m_factorizationIsOk = false; 00643 return; 00644 } 00645 00646 // Copy the U-segments to ucol(*) 00647 info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 00648 if ( info ) 00649 { 00650 m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() "; 00651 m_info = NumericalIssue; 00652 m_factorizationIsOk = false; 00653 return; 00654 } 00655 00656 // Form the L-segment 00657 info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); 00658 if ( info ) 00659 { 00660 m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT "; 00661 std::ostringstream returnInfo; 00662 returnInfo << info; 00663 m_lastError += returnInfo.str(); 00664 m_info = NumericalIssue; 00665 m_factorizationIsOk = false; 00666 return; 00667 } 00668 00669 // Update the determinant of the row permutation matrix 00670 // FIXME: the following test is not correct, we should probably take iperm_c into account and pivrow is not directly the row pivot. 00671 if (pivrow != jj) m_detPermR = -m_detPermR; 00672 00673 // Prune columns (0:jj-1) using column jj 00674 Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 00675 00676 // Reset repfnz for this column 00677 for (i = 0; i < nseg; i++) 00678 { 00679 irep = segrep(i); 00680 repfnz_k(irep) = emptyIdxLU; 00681 } 00682 } // end SparseLU within the panel 00683 jcol += panel_size; // Move to the next panel 00684 } // end for -- end elimination 00685 00686 m_detPermR = m_perm_r.determinant(); 00687 m_detPermC = m_perm_c.determinant(); 00688 00689 // Count the number of nonzeros in factors 00690 Base::countnz(n, m_nnzL, m_nnzU, m_glu); 00691 // Apply permutation to the L subscripts 00692 Base::fixupL(n, m_perm_r.indices(), m_glu); 00693 00694 // Create supernode matrix L 00695 m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); 00696 // Create the column major upper sparse matrix U; 00697 new (&m_Ustore) MappedSparseMatrix<Scalar, ColMajor, Index> ( m, n, m_nnzU, m_glu.xusub.data(), m_glu.usub.data(), m_glu.ucol.data() ); 00698 00699 m_info = Success; 00700 m_factorizationIsOk = true; 00701 } 00702 00703 template<typename MappedSupernodalType> 00704 struct SparseLUMatrixLReturnType : internal::no_assignment_operator 00705 { 00706 typedef typename MappedSupernodalType::Index Index; 00707 typedef typename MappedSupernodalType::Scalar Scalar; 00708 SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL) 00709 { } 00710 Index rows() { return m_mapL.rows(); } 00711 Index cols() { return m_mapL.cols(); } 00712 template<typename Dest> 00713 void solveInPlace( MatrixBase<Dest> &X) const 00714 { 00715 m_mapL.solveInPlace(X); 00716 } 00717 const MappedSupernodalType& m_mapL; 00718 }; 00719 00720 template<typename MatrixLType, typename MatrixUType> 00721 struct SparseLUMatrixUReturnType : internal::no_assignment_operator 00722 { 00723 typedef typename MatrixLType::Index Index; 00724 typedef typename MatrixLType::Scalar Scalar; 00725 SparseLUMatrixUReturnType(const MatrixLType& mapL, const MatrixUType& mapU) 00726 : m_mapL(mapL),m_mapU(mapU) 00727 { } 00728 Index rows() { return m_mapL.rows(); } 00729 Index cols() { return m_mapL.cols(); } 00730 00731 template<typename Dest> void solveInPlace(MatrixBase<Dest> &X) const 00732 { 00733 Index nrhs = X.cols(); 00734 Index n = X.rows(); 00735 // Backward solve with U 00736 for (Index k = m_mapL.nsuper(); k >= 0; k--) 00737 { 00738 Index fsupc = m_mapL.supToCol()[k]; 00739 Index lda = m_mapL.colIndexPtr()[fsupc+1] - m_mapL.colIndexPtr()[fsupc]; // leading dimension 00740 Index nsupc = m_mapL.supToCol()[k+1] - fsupc; 00741 Index luptr = m_mapL.colIndexPtr()[fsupc]; 00742 00743 if (nsupc == 1) 00744 { 00745 for (Index j = 0; j < nrhs; j++) 00746 { 00747 X(fsupc, j) /= m_mapL.valuePtr()[luptr]; 00748 } 00749 } 00750 else 00751 { 00752 Map<const Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(m_mapL.valuePtr()[luptr]), nsupc, nsupc, OuterStride<>(lda) ); 00753 Map< Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > U (&(X(fsupc,0)), nsupc, nrhs, OuterStride<>(n) ); 00754 U = A.template triangularView<Upper>().solve(U); 00755 } 00756 00757 for (Index j = 0; j < nrhs; ++j) 00758 { 00759 for (Index jcol = fsupc; jcol < fsupc + nsupc; jcol++) 00760 { 00761 typename MatrixUType::InnerIterator it(m_mapU, jcol); 00762 for ( ; it; ++it) 00763 { 00764 Index irow = it.index(); 00765 X(irow, j) -= X(jcol, j) * it.value(); 00766 } 00767 } 00768 } 00769 } // End For U-solve 00770 } 00771 const MatrixLType& m_mapL; 00772 const MatrixUType& m_mapU; 00773 }; 00774 00775 namespace internal { 00776 00777 template<typename _MatrixType, typename Derived, typename Rhs> 00778 struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 00779 : solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 00780 { 00781 typedef SparseLU<_MatrixType,Derived> Dec; 00782 EIGEN_MAKE_SOLVE_HELPERS(Dec,Rhs) 00783 00784 template<typename Dest> void evalTo(Dest& dst) const 00785 { 00786 dec()._solve(rhs(),dst); 00787 } 00788 }; 00789 00790 template<typename _MatrixType, typename Derived, typename Rhs> 00791 struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs> 00792 : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs> 00793 { 00794 typedef SparseLU<_MatrixType,Derived> Dec; 00795 EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs) 00796 00797 template<typename Dest> void evalTo(Dest& dst) const 00798 { 00799 this->defaultEvalTo(dst); 00800 } 00801 }; 00802 } // end namespace internal 00803 00804 } // End namespace Eigen 00805 00806 #endif