SparseLU.h
Go to the documentation of this file.
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) 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       // 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       for(Index j = 0; j < B.cols(); ++j)
00233         X.col(j) = rowsPermutation() * B.col(j);
00234       
00235       //Forward substitution with L
00236       this->matrixL().solveInPlace(X);
00237       this->matrixU().solveInPlace(X);
00238       
00239       // Permute back the solution 
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       // Initialize with the determinant of the row matrix
00260       Scalar det = Scalar(1.);
00261       //Note that the diagonal blocks of U are stored in supernodes,
00262       // which are available in the  L part :)
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     // Functions 
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     // Variables 
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; // The input (permuted ) matrix 
00334     SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
00335     MappedSparseMatrix<Scalar,ColMajor,Index> m_Ustore; // The upper triangular matrix
00336     PermutationType m_perm_c; // Column permutation 
00337     PermutationType m_perm_r ; // Row permutation
00338     IndexVector m_etree; // Column elimination tree 
00339     
00340     typename Base::GlobalLU_t m_glu; 
00341                                
00342     // SparseLU options 
00343     bool m_symmetricmode;
00344     // values for performance 
00345     internal::perfvalues<Index> m_perfv; 
00346     RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
00347     Index m_nnzL, m_nnzU; // Nonzeros in L and U factors 
00348     Index m_detPermR; // Determinant of the coefficient matrix
00349   private:
00350     // Disable copy constructor 
00351     SparseLU (const SparseLU& );
00352   
00353 }; // End class SparseLU
00354 
00355 
00356 
00357 // Functions needed by the anaysis phase
00368 template <typename MatrixType, typename OrderingType>
00369 void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
00370 {
00371   
00372   //TODO  It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
00373   
00374   OrderingType ord; 
00375   ord(mat,m_perm_c);
00376   
00377   // Apply the permutation to the column of the input  matrix
00378   //First copy the whole input matrix. 
00379   m_mat = mat;
00380   if (m_perm_c.size()) {
00381     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.  
00382     //Then, permute only the column pointers
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   // Compute the column elimination tree of the permuted matrix 
00399   IndexVector firstRowElt;
00400   internal::coletree(m_mat, m_etree,firstRowElt); 
00401      
00402   // In symmetric mode, do not do postorder here
00403   if (!m_symmetricmode) {
00404     IndexVector post, iwork; 
00405     // Post order etree
00406     internal::treePostorder(m_mat.cols(), m_etree, post); 
00407       
00408    
00409     // Renumber etree in postorder 
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     // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
00416     PermutationType post_perm(m); 
00417     for (Index i = 0; i < m; i++) 
00418       post_perm.indices()(i) = post(i); 
00419         
00420     // Combine the two permutations : postorder the permutation for future use
00421     if(m_perm_c.size()) {
00422       m_perm_c = post_perm * m_perm_c;
00423     }
00424     
00425   } // end postordering 
00426   
00427   m_analysisIsOk = true; 
00428 }
00429 
00430 // Functions needed by the numerical factorization phase
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   // Apply the column permutation computed in analyzepattern()
00462   //   m_mat = matrix * m_perm_c.inverse(); 
00463   m_mat = matrix;
00464   if (m_perm_c.size()) 
00465   {
00466     m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
00467     //Then, permute only the column pointers
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   { //FIXME This should not be needed if the empty permutation is handled transparently
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   // Allocate working storage common to the factor routines
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   // Set up pointers for integer working arrays 
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   // Set up pointers for scalar working arrays 
00516   ScalarVector dense; 
00517   dense.setZero(maxpanel);
00518   ScalarVector tempv; 
00519   tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
00520   
00521   // Compute the inverse of perm_c
00522   PermutationType iperm_c(m_perm_c.inverse()); 
00523   
00524   // Identify initial relaxed snodes
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; // Record the determinant of the row permutation
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   // Work on one 'panel' at a time. A panel is one of the following :
00541   //  (a) a relaxed supernode at the bottom of the etree, or
00542   //  (b) panel_size contiguous columns, <panel_size> defined by the user
00543   Index jcol; 
00544   IndexVector panel_histo(n);
00545   Index pivrow; // Pivotal row number in the original row matrix
00546   Index nseg1; // Number of segments in U-column above panel row jcol
00547   Index nseg; // Number of segments in each U-column 
00548   Index irep; 
00549   Index i, k, jj; 
00550   for (jcol = 0; jcol < n; )
00551   {
00552     // Adjust panel size so that a panel won't overlap with the next relaxed snode. 
00553     Index panel_size = m_perfv.panel_size; // upper bound on panel width
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     // Symbolic outer factorization on a panel of columns 
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     // Numeric sup-panel updates in topological order 
00569     Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 
00570     
00571     // Sparse LU within the panel, and below the panel diagonal 
00572     for ( jj = jcol; jj< jcol + panel_size; jj++) 
00573     {
00574       k = (jj - jcol) * m; // Column index for w-wide arrays 
00575       
00576       nseg = nseg1; // begin after all the panel segments
00577       //Depth-first-search for the current column
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       // Numeric updates to this column 
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       // Copy the U-segments to ucol(*)
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       // Form the L-segment 
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       // Update the determinant of the row permutation matrix
00624       if (pivrow != jj) m_detPermR *= -1;
00625 
00626       // Prune columns (0:jj-1) using column jj
00627       Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
00628       
00629       // Reset repfnz for this column 
00630       for (i = 0; i < nseg; i++)
00631       {
00632         irep = segrep(i); 
00633         repfnz_k(irep) = emptyIdxLU; 
00634       }
00635     } // end SparseLU within the panel  
00636     jcol += panel_size;  // Move to the next panel
00637   } // end for -- end elimination 
00638   
00639   // Count the number of nonzeros in factors 
00640   Base::countnz(n, m_nnzL, m_nnzU, m_glu); 
00641   // Apply permutation  to the L subscripts 
00642   Base::fixupL(n, m_perm_r.indices(), m_glu); 
00643   
00644   // Create supernode matrix L 
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   // Create the column major upper sparse matrix  U; 
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     // Backward solve with U
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]; // leading dimension
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     } // End For U-solve
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 } // end namespace internal
00753 
00754 } // End namespace Eigen 
00755 
00756 #endif


acado
Author(s): Milan Vukov, Rien Quirynen
autogenerated on Thu Aug 27 2015 12:00:33