Add additional methods in SparseLU and Improve the naming conventions

1 files changed, 126 insertions, 74 deletions

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 175794811..b1a74581a 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -2,6 +2,7 @@
 // for linear algebra.
 //
 // Copyright (C) 2012 Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
+// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@inria.fr>
 //
 // This Source Code Form is subject to the terms of the Mozilla
 // Public License v. 2.0. If a copy of the MPL was not distributed
@@ -13,6 +14,8 @@
 
 namespace Eigen {
 
+template <typename _MatrixType, typename _OrderingType> class SparseLU;
+template <typename MappedSparseMatrixType> struct SparseLUMatrixLReturnType; 
 /** \ingroup SparseLU_Module
   * \class SparseLU
   * 
@@ -65,7 +68,7 @@ namespace Eigen {
   * \sa \ref OrderingMethods_Module
   */
 template <typename _MatrixType, typename _OrderingType>
-class SparseLU
+class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typename _MatrixType::Index>
 {
   public:
     typedef _MatrixType MatrixType; 
@@ -74,17 +77,18 @@ class SparseLU
     typedef typename MatrixType::RealScalar RealScalar; 
     typedef typename MatrixType::Index Index; 
     typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix;
-    typedef SuperNodalMatrix<Scalar, Index> SCMatrix; 
+    typedef internal::MappedSuperNodalMatrix<Scalar, Index> SCMatrix; 
     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
     typedef Matrix<Index,Dynamic,1> IndexVector;
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+    typedef internal::SparseLUImpl<Scalar, Index> Base;
     
   public:
-    SparseLU():m_isInitialized(true),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
+    SparseLU():m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
     {
       initperfvalues(); 
     }
-    SparseLU(const MatrixType& matrix):m_isInitialized(true),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
+    SparseLU(const MatrixType& matrix):m_isInitialized(true),m_lastError(""),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
     {
       initperfvalues(); 
       compute(matrix);
@@ -119,34 +123,22 @@ class SparseLU
       m_symmetricmode = sym;
     }
     
-    /** Set the threshold used for a diagonal entry to be an acceptable pivot. */
-    void diagPivotThresh(RealScalar thresh)
+    /** Returns an expression of the matrix L, internally stored as supernodes 
+     * For a triangular solve with this matrix, use
+     * \code
+     * y = b; matrixL().solveInPlace(y);
+     * \endcode
+     */
+    SparseLUMatrixLReturnType<SCMatrix> matrixL() const
     {
-      m_diagpivotthresh = thresh; 
+      return SparseLUMatrixLReturnType<SCMatrix>(m_Lstore);
     }
-     
-    /** Return the number of nonzero elements in the L factor */
-    int nnzL()
-    {
-      if (m_factorizationIsOk)
-        return m_nnzL; 
-      else
-      {
-        std::cerr<<"Numerical factorization should be done before\n"; 
-        return 0; 
-      }
-    }
-    /** Return the number of nonzero elements in the U factor */
-    int nnzU()
+    /** Set the threshold used for a diagonal entry to be an acceptable pivot. */
+    void setPivotThreshold(RealScalar thresh)
     {
-      if (m_factorizationIsOk)
-        return m_nnzU; 
-      else
-      {
-        std::cerr<<"Numerical factorization should be done before\n"; 
-        return 0; 
-      }
+      m_diagpivotthresh = thresh; 
     }
+
     /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
       *
       * \sa compute()
@@ -160,6 +152,18 @@ class SparseLU
           return internal::solve_retval<SparseLU, Rhs>(*this, B.derived());
     }
 
+        /** \returns the solution X of \f$ A X = B \f$ using the current decomposition of A.
+      *
+      * \sa compute()
+      */
+    template<typename Rhs>
+    inline const internal::sparse_solve_retval<SparseLU, Rhs> solve(const SparseMatrixBase<Rhs>& B) const 
+    {
+      eigen_assert(m_factorizationIsOk && "SparseLU is not initialized."); 
+      eigen_assert(rows()==B.rows()
+                    && "SparseLU::solve(): invalid number of rows of the right hand side matrix B");
+          return internal::sparse_solve_retval<SparseLU, Rhs>(*this, B.derived());
+    }
     
      /** \brief Reports whether previous computation was successful.
       *
@@ -174,7 +178,13 @@ class SparseLU
       eigen_assert(m_isInitialized && "Decomposition is not initialized.");
       return m_info;
     }
-
+    /**
+     * \returns A string describing the type of error
+     */
+    std::string lastErrorMessage() const
+    {
+      return m_lastError; 
+    }
     template<typename Rhs, typename Dest>
     bool _solve(const MatrixBase<Rhs> &B, MatrixBase<Dest> &_X) const
     {
@@ -194,7 +204,8 @@ class SparseLU
         X.col(j) = m_perm_r * B.col(j); 
       
       //Forward substitution with L 
-      m_Lstore.solveInPlace(X);
+//       m_Lstore.solveInPlace(X);
+        this->matrixL().solveInPlace(X);
       
       // Backward solve with U
       for (int k = m_Lstore.nsuper(); k >= 0; k--)
@@ -256,6 +267,7 @@ class SparseLU
     bool m_isInitialized;
     bool m_factorizationIsOk;
     bool m_analysisIsOk;
+    std::string m_lastError;
     NCMatrix m_mat; // The input (permuted ) matrix 
     SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
     MappedSparseMatrix<Scalar> m_Ustore; // The upper triangular matrix
@@ -263,16 +275,15 @@ class SparseLU
     PermutationType m_perm_r ; // Row permutation
     IndexVector m_etree; // Column elimination tree 
     
-    LU_GlobalLU_t<IndexVector, ScalarVector> m_glu; 
+    typename Base::GlobalLU_t m_glu; 
                                
-    // SuperLU/SparseLU options 
+    // SparseLU options 
     bool m_symmetricmode;
-    
     // values for performance 
-    LU_perfvalues m_perfv; 
+    internal::perfvalues m_perfv; 
     RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
     int m_nnzL, m_nnzU; // Nonzeros in L and U factors 
-  
+    
   private:
     // Copy constructor 
     SparseLU (SparseLU& ) {}
@@ -301,18 +312,17 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   ord(mat,m_perm_c);
   
   // Apply the permutation to the column of the input  matrix
-//   m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix
-   
   //First copy the whole input matrix. 
   m_mat = mat;
-  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.  
-  //Then, permute only the column pointers
-  for (int i = 0; i < mat.cols(); i++)
-  {
-    m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i]; 
-    m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 
+  if (m_perm_c.size()) {
+    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.  
+    //Then, permute only the column pointers
+    for (int i = 0; i < mat.cols(); i++)
+    {
+      m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i]; 
+      m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 
+    }
   }
-    
   // Compute the column elimination tree of the permuted matrix 
   IndexVector firstRowElt;
   internal::coletree(m_mat, m_etree,firstRowElt); 
@@ -331,12 +341,14 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
     m_etree = iwork;
     
     // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
-    PermutationType post_perm(m); //FIXME Use directly a constructor with post
+    PermutationType post_perm(m); 
     for (int i = 0; i < m; i++) 
       post_perm.indices()(i) = post(i); 
         
     // Combine the two permutations : postorder the permutation for future use
-    m_perm_c = post_perm * m_perm_c;
+    if(m_perm_c.size()) {
+      m_perm_c = post_perm * m_perm_c;
+    }
     
   } // end postordering 
   
@@ -367,7 +379,7 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
 template <typename MatrixType, typename OrderingType>
 void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
 {
-  
+  using internal::emptyIdxLU;
   eigen_assert(m_analysisIsOk && "analyzePattern() should be called first"); 
   eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
   
@@ -377,12 +389,20 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   // Apply the column permutation computed in analyzepattern()
   //   m_mat = matrix * m_perm_c.inverse(); 
   m_mat = matrix;
-  m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
-  //Then, permute only the column pointers
-  for (int i = 0; i < matrix.cols(); i++)
+  if (m_perm_c.size()) 
   {
-    m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i]; 
-    m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i+1] - matrix.outerIndexPtr()[i]; 
+    m_mat.uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
+    //Then, permute only the column pointers
+    for (int i = 0; i < matrix.cols(); i++)
+    {
+      m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i]; 
+      m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i+1] - matrix.outerIndexPtr()[i]; 
+    }
+  } 
+  else 
+  { //FIXME This should not be needed if the empty permutation is handled transparently
+    m_perm_c.resize(matrix.cols());
+    for(int i = 0; i < matrix.cols(); ++i) m_perm_c.indices()(i) = i;
   }
   
   int m = m_mat.rows();
@@ -391,10 +411,10 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   int maxpanel = m_perfv.panel_size * m;
   // Allocate working storage common to the factor routines
   int lwork = 0;
-  int info = SparseLUBase<Scalar,Index>::LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 
+  int info = Base::memInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 
   if (info) 
   {
-    std::cerr << "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
+    m_lastError = "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
     m_factorizationIsOk = false;
     return ; 
   }
@@ -406,7 +426,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   IndexVector repfnz(maxpanel);
   IndexVector panel_lsub(maxpanel);
   IndexVector xprune(n); xprune.setZero();
-  IndexVector marker(m*LU_NO_MARKER); marker.setZero();
+  IndexVector marker(m*internal::LUNoMarker); marker.setZero();
   
   repfnz.setConstant(-1); 
   panel_lsub.setConstant(-1);
@@ -415,7 +435,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   ScalarVector dense; 
   dense.setZero(maxpanel);
   ScalarVector tempv; 
-  tempv.setZero(LU_NUM_TEMPV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
+  tempv.setZero(internal::LUnumTempV(m, m_perfv.panel_size, m_perfv.maxsuper, /*m_perfv.rowblk*/m) );
   
   // Compute the inverse of perm_c
   PermutationType iperm_c(m_perm_c.inverse()); 
@@ -423,16 +443,16 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   // Identify initial relaxed snodes
   IndexVector relax_end(n);
   if ( m_symmetricmode == true ) 
-    SparseLUBase<Scalar,Index>::LU_heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
+    Base::heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
   else
-    SparseLUBase<Scalar,Index>::LU_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
+    Base::relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
   
   
   m_perm_r.resize(m); 
   m_perm_r.indices().setConstant(-1);
   marker.setConstant(-1);
   
-  m_glu.supno(0) = IND_EMPTY; m_glu.xsup.setConstant(0);
+  m_glu.supno(0) = emptyIdxLU; m_glu.xsup.setConstant(0);
   m_glu.xsup(0) = m_glu.xlsub(0) = m_glu.xusub(0) = m_glu.xlusup(0) = Index(0);
   
   // Work on one 'panel' at a time. A panel is one of the following :
@@ -451,7 +471,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
     int panel_size = m_perfv.panel_size; // upper bound on panel width
     for (k = jcol + 1; k < (std::min)(jcol+panel_size, n); k++)
     {
-      if (relax_end(k) != IND_EMPTY) 
+      if (relax_end(k) != emptyIdxLU) 
       {
         panel_size = k - jcol; 
         break; 
@@ -461,10 +481,10 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       panel_size = n - jcol; 
       
     // Symbolic outer factorization on a panel of columns 
-    SparseLUBase<Scalar,Index>::LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r.indices(), nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_glu); 
+    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); 
     
     // Numeric sup-panel updates in topological order 
-    SparseLUBase<Scalar,Index>::LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 
+    Base::panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_glu); 
     
     // Sparse LU within the panel, and below the panel diagonal 
     for ( jj = jcol; jj< jcol + panel_size; jj++) 
@@ -475,10 +495,10 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       //Depth-first-search for the current column
       VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
       VectorBlock<IndexVector> repfnz_k(repfnz, k, m); 
-      info = SparseLUBase<Scalar,Index>::LU_column_dfs(m, jj, m_perm_r.indices(), m_perfv.maxsuper, nseg, panel_lsubk, segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 
+      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); 
       if ( info ) 
       {
-        std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n";
+        m_lastError =  "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
@@ -486,52 +506,55 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       // Numeric updates to this column 
       VectorBlock<ScalarVector> dense_k(dense, k, m); 
       VectorBlock<IndexVector> segrep_k(segrep, nseg1, m-nseg1); 
-      info = SparseLUBase<Scalar,Index>::LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 
+      info = Base::column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 
       if ( info ) 
       {
-        std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() \n";
+        m_lastError = "UNABLE TO EXPAND MEMORY IN COLUMN_BMOD() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Copy the U-segments to ucol(*)
-      info = SparseLUBase<Scalar,Index>::LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 
+      info = Base::copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 
       if ( info ) 
       {
-        std::cerr << "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() \n";
+        m_lastError = "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() ";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Form the L-segment 
-      info = SparseLUBase<Scalar,Index>::LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
+      info = Base::pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
       if ( info ) 
       {
-        std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl; 
+        m_lastError = "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT ";
+        std::ostringstream returnInfo;
+        returnInfo << info; 
+        m_lastError += returnInfo.str();
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
         return; 
       }
       
       // Prune columns (0:jj-1) using column jj
-      SparseLUBase<Scalar,Index>::LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
+      Base::pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
       
       // Reset repfnz for this column 
       for (i = 0; i < nseg; i++)
       {
         irep = segrep(i); 
-        repfnz_k(irep) = IND_EMPTY; 
+        repfnz_k(irep) = emptyIdxLU; 
       }
     } // end SparseLU within the panel  
     jcol += panel_size;  // Move to the next panel
   } // end for -- end elimination 
   
   // Count the number of nonzeros in factors 
-  SparseLUBase<Scalar,Index>::LU_countnz(n, m_nnzL, m_nnzU, m_glu); 
+  Base::countnz(n, m_nnzL, m_nnzU, m_glu); 
   // Apply permutation  to the L subscripts 
-  SparseLUBase<Scalar,Index>::LU_fixupL(n, m_perm_r.indices(), m_glu); 
+  Base::fixupL(n, m_perm_r.indices(), m_glu); 
   
   // Create supernode matrix L 
   m_Lstore.setInfos(m, n, m_glu.lusup, m_glu.xlusup, m_glu.lsub, m_glu.xlsub, m_glu.supno, m_glu.xsup); 
@@ -542,6 +565,23 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   m_factorizationIsOk = true;
 }
 
+template<typename MappedSupernodalType>
+struct SparseLUMatrixLReturnType
+{
+  typedef typename MappedSupernodalType::Index Index;
+  typedef typename MappedSupernodalType::Scalar Scalar;
+  SparseLUMatrixLReturnType(const MappedSupernodalType& mapL) : m_mapL(mapL)
+  { }
+  Index rows() { return m_mapL.rows(); }
+  Index cols() { return m_mapL.cols(); }
+  template<typename Dest>
+  void solveInPlace( MatrixBase<Dest> &X) const
+  {
+    m_mapL.solveInPlace(X);
+  }
+  const MappedSupernodalType& m_mapL;
+};
+
 namespace internal {
   
 template<typename _MatrixType, typename Derived, typename Rhs>
@@ -557,6 +597,18 @@ struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
   }
 };
 
+template<typename _MatrixType, typename Derived, typename Rhs>
+struct sparse_solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
+  : sparse_solve_retval_base<SparseLU<_MatrixType,Derived>, Rhs>
+{
+  typedef SparseLU<_MatrixType,Derived> Dec;
+  EIGEN_MAKE_SPARSE_SOLVE_HELPERS(Dec,Rhs)
+
+  template<typename Dest> void evalTo(Dest& dst) const
+  {
+    this->defaultEvalTo(dst);
+  }
+};
 } // end namespace internal
 
 } // End namespace Eigen