Define sparseLU functions as static

16 files changed, 130 insertions, 179 deletions

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 77df091c3..f5d15ec6b 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -18,6 +18,8 @@ namespace Eigen {
 #include "SparseLU_Structs.h"
 #include "SparseLU_Matrix.h"
 
+// Base structure containing all the factorization routines
+#include "SparseLUBase.h"
 /**
  * \ingroup SparseLU_Module
  * \brief Sparse supernodal LU factorization for general matrices
@@ -40,6 +42,7 @@ class SparseLU
     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
     typedef Matrix<Index,Dynamic,1> IndexVector;
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+    
   public:
     SparseLU():m_isInitialized(true),m_Ustore(0,0,0,0,0,0),m_symmetricmode(false),m_diagpivotthresh(1.0)
     {
@@ -58,6 +61,7 @@ class SparseLU
     
     void analyzePattern (const MatrixType& matrix);
     void factorize (const MatrixType& matrix);
+    void simplicialfactorize(const MatrixType& matrix);
     
     /**
      * Compute the symbolic and numeric factorization of the input sparse matrix.
@@ -224,8 +228,7 @@ class SparseLU
     PermutationType m_perm_r ; // Row permutation
     IndexVector m_etree; // Column elimination tree 
     
-    LU_GlobalLU_t<IndexVector, ScalarVector> m_glu; // persistent data to facilitate multiple factors 
-                               // FIXME All fields of this struct can be defined separately as class members
+    LU_GlobalLU_t<IndexVector, ScalarVector> m_glu; 
                                
     // SuperLU/SparseLU options 
     bool m_symmetricmode;
@@ -243,7 +246,6 @@ class SparseLU
 
 
 // Functions needed by the anaysis phase
-#include "SparseLU_Coletree.h"
 /** 
  * Compute the column permutation to minimize the fill-in (file amd.c )
  * 
@@ -262,9 +264,6 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   
   OrderingType ord; 
   ord(mat,m_perm_c);
-  //FIXME Check the right semantic behind m_perm_c
-  // that is, column j of mat goes to column m_perm_c(j) of 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
@@ -282,13 +281,13 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   // Compute the column elimination tree of the permuted matrix 
   /*if (m_etree.size() == 0)  */m_etree.resize(m_mat.cols());
   
-  LU_sp_coletree(m_mat, m_etree); 
+  SparseLUBase<Scalar,Index>::LU_sp_coletree(m_mat, m_etree); 
      
   // In symmetric mode, do not do postorder here
   if (!m_symmetricmode) {
     IndexVector post, iwork; 
     // Post order etree
-    LU_TreePostorder(m_mat.cols(), m_etree, post); 
+    SparseLUBase<Scalar,Index>::LU_TreePostorder(m_mat.cols(), m_etree, post); 
       
    
     // Renumber etree in postorder 
@@ -310,21 +309,7 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   m_analysisIsOk = true; 
 }
 
-// Functions needed by the numerical factorization phase 
-#include "SparseLU_Memory.h"
-#include "SparseLU_heap_relax_snode.h"
-#include "SparseLU_relax_snode.h"
-#include "SparseLU_snode_dfs.h"
-#include "SparseLU_snode_bmod.h"
-#include "SparseLU_pivotL.h"
-#include "SparseLU_panel_dfs.h"
-#include "SparseLU_kernel_bmod.h"
-#include "SparseLU_panel_bmod.h"
-#include "SparseLU_column_dfs.h"
-#include "SparseLU_column_bmod.h"
-#include "SparseLU_copy_to_ucol.h"
-#include "SparseLU_pruneL.h"
-#include "SparseLU_Utils.h"
+// Functions needed by the numerical factorization phase
 
 
 /** 
@@ -370,7 +355,7 @@ 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 = LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 
+  int info = SparseLUBase<Scalar,Index>::LUMemInit(m, n, nnz, lwork, m_perfv.fillfactor, m_perfv.panel_size, m_glu); 
   if (info) 
   {
     std::cerr << "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
@@ -402,25 +387,17 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   // Identify initial relaxed snodes
   IndexVector relax_end(n);
   if ( m_symmetricmode == true ) 
-    LU_heap_relax_snode<IndexVector>(n, m_etree, m_perfv.relax, marker, relax_end);
+    SparseLUBase<Scalar,Index>::LU_heap_relax_snode(n, m_etree, m_perfv.relax, marker, relax_end);
   else
-    LU_relax_snode<IndexVector>(n, m_etree, m_perfv.relax, marker, relax_end);
+    SparseLUBase<Scalar,Index>::LU_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);
   
-  IndexVector& xsup = m_glu.xsup; 
-  IndexVector& supno = m_glu.supno; 
-  IndexVector& xlsub = m_glu.xlsub;
-  IndexVector& xlusup = m_glu.xlusup;
-  IndexVector& xusub = m_glu.xusub;
-  ScalarVector& lusup = m_glu.lusup; 
-  Index& nzlumax = m_glu.nzlumax; 
-    
-  supno(0) = IND_EMPTY; xsup.setConstant(0);
-  xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = Index(0);
+  m_glu.supno(0) = IND_EMPTY; 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 :
   //  (a) a relaxed supernode at the bottom of the etree, or
@@ -441,7 +418,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       
       // Factorize the relaxed supernode(jcol:kcol)
       // First, determine the union of the row structure of the snode 
-      info = LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); 
+      info = SparseLUBase<Scalar,Index>::LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); 
       if ( info ) 
       {
         std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n";
@@ -449,15 +426,15 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
         m_factorizationIsOk = false; 
         return; 
       }
-      nextu = xusub(jcol); //starting location of column jcol in ucol
-      nextlu = xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes)
-      jsupno = supno(jcol); // Supernode number which column jcol belongs to 
-      fsupc = xsup(jsupno); //First column number of the current supernode
-      new_next = nextlu + (xlsub(fsupc+1)-xlsub(fsupc)) * (kcol - jcol + 1);
+      nextu = m_glu.xusub(jcol); //starting location of column jcol in ucol
+      nextlu = m_glu.xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes)
+      jsupno = m_glu.supno(jcol); // Supernode number which column jcol belongs to 
+      fsupc = m_glu.xsup(jsupno); //First column number of the current supernode
+      new_next = nextlu + (m_glu.xlsub(fsupc+1)-m_glu.xlsub(fsupc)) * (kcol - jcol + 1);
       int mem; 
-      while (new_next > nzlumax ) 
+      while (new_next > m_glu.nzlumax ) 
       {
-        mem = LUMemXpand(lusup, nzlumax, nextlu, LUSUP, m_glu.num_expansions);
+        mem = SparseLUBase<Scalar,Index>::LUMemXpand(m_glu.lusup, m_glu.nzlumax, nextlu, LUSUP, m_glu.num_expansions);
         if (mem) 
         {
           std::cerr << "MEMORY ALLOCATION FAILED FOR L FACTOR \n"; 
@@ -468,16 +445,16 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       
       // Now, left-looking factorize each column within the snode
       for (icol = jcol; icol<=kcol; icol++){
-        xusub(icol+1) = nextu;
+        m_glu.xusub(icol+1) = nextu;
         // Scatter into SPA dense(*)
         for (typename MatrixType::InnerIterator it(m_mat, icol); it; ++it)
           dense(it.row()) = it.value();
         
         // Numeric update within the snode 
-        LU_snode_bmod(icol, fsupc, dense, m_glu); 
+        SparseLUBase<Scalar,Index>::LU_snode_bmod(icol, fsupc, dense, m_glu); 
         
         // Eliminate the current column 
-        info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); 
+        info = SparseLUBase<Scalar,Index>::LU_pivotL(icol, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu); 
         if ( info ) 
         {
           m_info = NumericalIssue; 
@@ -505,10 +482,10 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
         panel_size = n - jcol; 
         
       // Symbolic outer factorization on a panel of columns 
-      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); 
+      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); 
       
       // Numeric sup-panel updates in topological order 
-      LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu); 
+      SparseLUBase<Scalar,Index>::LU_panel_bmod(m, panel_size, jcol, nseg1, dense, tempv, segrep, repfnz, m_perfv, m_glu); 
       
       // Sparse LU within the panel, and below the panel diagonal 
       for ( jj = jcol; jj< jcol + panel_size; jj++) 
@@ -519,7 +496,7 @@ 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 = 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 = 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); 
         if ( info ) 
         {
           std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n";
@@ -530,7 +507,7 @@ 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 = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 
+        info = SparseLUBase<Scalar,Index>::LU_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";
@@ -540,7 +517,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
         }
         
         // Copy the U-segments to ucol(*)
-        info = LU_copy_to_ucol(jj, nseg, segrep, repfnz_k ,m_perm_r.indices(), dense_k, m_glu); 
+        info = SparseLUBase<Scalar,Index>::LU_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";
@@ -550,7 +527,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
         }
         
         // Form the L-segment 
-        info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r.indices(), iperm_c.indices(), pivrow, m_glu);
+        info = SparseLUBase<Scalar,Index>::LU_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; 
@@ -560,7 +537,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
         }
         
         // Prune columns (0:jj-1) using column jj
-        LU_pruneL(jj, m_perm_r.indices(), pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
+        SparseLUBase<Scalar,Index>::LU_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++)
@@ -574,11 +551,9 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   } // end for -- end elimination 
   
   // Count the number of nonzeros in factors 
-  LU_countnz(n, m_nnzL, m_nnzU, m_glu); 
+  SparseLUBase<Scalar,Index>::LU_countnz(n, m_nnzL, m_nnzU, m_glu); 
   // Apply permutation  to the L subscripts 
-  LU_fixupL(n, m_perm_r.indices(), m_glu); 
-  
-  
+  SparseLUBase<Scalar,Index>::LU_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); 
@@ -589,7 +564,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   m_factorizationIsOk = true;
 }
 
-
+// #include "SparseLU_simplicialfactorize.h"
 namespace internal {
   
 template<typename _MatrixType, typename Derived, typename Rhs>
@@ -607,7 +582,5 @@ struct solve_retval<SparseLU<_MatrixType,Derived>, Rhs>
 
 } // end namespace internal
 
-
-
 } // End namespace Eigen 
-#endif
\ No newline at end of file
+#endif
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
index 964f5e433..bb4067a45 100644
--- a/Eigen/src/SparseLU/SparseLU_Coletree.h
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -31,8 +31,8 @@
 #ifndef SPARSELU_COLETREE_H
 #define SPARSELU_COLETREE_H
 /** Find the root of the tree/set containing the vertex i : Use Path halving */ 
-template<typename IndexVector>
-int etree_find (int i, IndexVector& pp)
+template< typename Scalar,typename Index>
+int SparseLUBase<Scalar,Index>::etree_find (int i, IndexVector& pp)
 {
   int p = pp(i); // Parent 
   int gp = pp(p); // Grand parent 
@@ -50,8 +50,8 @@ int etree_find (int i, IndexVector& pp)
   * NOTE : The matrix is supposed to be in column-major format. 
   * 
   */
-template<typename MatrixType, typename IndexVector>
-int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
 {
   int nc = mat.cols(); // Number of columns 
   int nr = mat.rows(); // Number of rows 
@@ -106,8 +106,8 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
   * Depth-first search from vertex n.  No recursion.
   * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
 */
-template<typename IndexVector>
-void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVector& next_kid, IndexVector& post, int postnum)
 {
   int current = n, first, next;
   while (postnum != n) 
@@ -152,8 +152,8 @@ void LU_nr_etdfs (int n, IndexVector& parent, IndexVector& first_kid, IndexVecto
   * \param parent Input tree
   * \param post postordered tree
   */
-template<typename IndexVector>
-void LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
 {
   IndexVector first_kid, next_kid; // Linked list of children 
   int postnum; 
diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h
index 48b36f5b4..0396ab61f 100644
--- a/Eigen/src/SparseLU/SparseLU_Memory.h
+++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -49,8 +49,9 @@
   * \param keep_prev  1: use length  and do not expand the vector; 0: compute new_len and expand
   * \param [in,out]num_expansions Number of times the memory has been expanded
   */
-template <typename VectorType >
-int  expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions) 
+template <typename Scalar, typename Index>
+template <typename VectorType>
+int  SparseLUBase<Scalar,Index>::expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_expansions) 
 {
   
   float alpha = 1.5; // Ratio of the memory increase 
@@ -122,12 +123,9 @@ int  expand(VectorType& vec, int& length, int nbElts, int keep_prev, int& num_ex
  * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated memory when allocation failed, and 0 on success
  * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the same row permutation
  */
-template <typename IndexVector,typename ScalarVector>
-int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,  LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,  GlobalLU_t& glu)
 {
-  typedef typename ScalarVector::Scalar Scalar; 
-  typedef typename IndexVector::Scalar Index; 
-  
   int& num_expansions = glu.num_expansions; //No memory expansions so far
   num_expansions = 0; 
   glu.nzumax = glu.nzlumax = std::max(fillratio * annz, m*n); // estimated number of nonzeros in U 
@@ -187,8 +185,9 @@ int LUMemInit(int m, int n, int annz, int lwork, int fillratio, int panel_size,
  * \param glu Global data structure 
  * \return 0 on success, > 0 size of the memory allocated so far
  */
+template <typename Scalar, typename Index>
 template <typename VectorType>
-int LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions)
+int SparseLUBase<Scalar,Index>::LUMemXpand(VectorType& vec, int& maxlen, int nbElts, LU_MemType memtype, int& num_expansions)
 {
   int failed_size; 
   if (memtype == USUB)
diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h
index 316b09ab0..b13930dbb 100644
--- a/Eigen/src/SparseLU/SparseLU_Utils.h
+++ b/Eigen/src/SparseLU/SparseLU_Utils.h
@@ -15,8 +15,8 @@
 /**
  * \brief Count Nonzero elements in the factors
  */
-template <typename IndexVector, typename ScalarVector>
-void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_countnz(const int n, int& nnzL, int& nnzU, GlobalLU_t& glu)
 {
  nnzL = 0; 
  nnzU = (glu.xusub)(n); 
@@ -46,8 +46,8 @@ void LU_countnz(const int n, int& nnzL, int& nnzU, LU_GlobalLU_t<IndexVector, Sc
  * and applies permutation to the remaining subscripts
  * 
  */
-template <typename IndexVector, typename ScalarVector>
-void LU_fixupL(const int n, const IndexVector& perm_r, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_fixupL(const int n, const IndexVector& perm_r, GlobalLU_t& glu)
 {
   int fsupc, i, j, k, jstart; 
   
diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h
index bf25a33fc..b268c4348 100644
--- a/Eigen/src/SparseLU/SparseLU_column_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
@@ -46,11 +46,9 @@
  *         > 0 - number of bytes allocated when run out of space
  * 
  */
-template <typename IndexVector, typename ScalarVector, typename BlockIndexVector, typename BlockScalarVector>
-int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, ScalarVector& tempv, BlockIndexVector& segrep, BlockIndexVector& repfnz, int fpanelc, GlobalLU_t& glu)
 {
-  typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar; 
   int  jsupno, k, ksub, krep, ksupno; 
   int lptr, nrow, isub, irow, nextlu, new_next, ufirst; 
   int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros; 
@@ -95,9 +93,6 @@ int LU_column_bmod(const int jcol, const int nseg, BlockScalarVector& dense, Sca
       nsupr = glu.xlsub(fsupc+1) - glu.xlsub(fsupc); 
       nrow = nsupr - d_fsupc - nsupc; 
       
-      // NOTE  Unlike the original implementation in SuperLU, the only feature  
-      // available here is a sup-col update. 
-      
       // Perform a triangular solver and block update, 
       // then scatter the result of sup-col update to dense
       no_zeros = kfnz - fst_col; 
diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h
index 568e0686c..fa8dcf18d 100644
--- a/Eigen/src/SparseLU/SparseLU_column_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h
@@ -42,15 +42,15 @@
  * \param m number of rows in the matrix
  * \param jcol Current column 
  * \param perm_r Row permutation
- * \param maxsuper 
+ * \param maxsuper  Maximum number of column allowed in a supernode
  * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
  * \param lsub_col defines the rhs vector to start the dfs
  * \param [in,out] segrep Segment representatives - new segments appended 
- * \param repfnz
+ * \param repfnz  First nonzero location in each row
  * \param xprune 
- * \param marker
+ * \param marker  marker[i] == jj, if i was visited during dfs of current column jj;
  * \param parent
- * \param xplore
+ * \param xplore working array
  * \param glu global LU data 
  * \return 0 success
  *         > 0 number of bytes allocated when run out of space
@@ -60,6 +60,7 @@ template<typename IndexVector, typename ScalarVector>
 struct LU_column_dfs_traits
 {
   typedef typename IndexVector::Scalar Index;
+  typedef typename ScalarVector::Scalar Scalar;
   LU_column_dfs_traits(Index jcol, Index& jsuper, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
    : m_jcol(jcol), m_jsuper_ref(jsuper), m_glu(glu)
  {}
@@ -70,7 +71,7 @@ struct LU_column_dfs_traits
   void mem_expand(IndexVector& lsub, int& nextl, int chmark)
   {
     if (nextl >= m_glu.nzlmax)
-      LUMemXpand<IndexVector>(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); 
+      SparseLUBase<Scalar,Index>::LUMemXpand(lsub, m_glu.nzlmax, nextl, LSUB, m_glu.num_expansions); 
     if (chmark != (m_jcol-1)) m_jsuper_ref = IND_EMPTY;
   }
   enum { ExpandMem = true };
@@ -80,11 +81,9 @@ struct LU_column_dfs_traits
   LU_GlobalLU_t<IndexVector, ScalarVector>& m_glu;
 };
 
-template <typename IndexVector, typename ScalarVector, typename BlockIndexVector>
-int LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg,  BlockIndexVector& lsub_col, IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, int maxsuper, int& nseg,  BlockIndexVector& lsub_col, IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
 {
-  typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar; 
   
   int jsuper = glu.supno(jcol); 
   int nextl = glu.xlsub(jcol); 
diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
index 541785881..d3227469d 100644
--- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
+++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
@@ -43,11 +43,9 @@
  *         > 0 - number of bytes allocated when run out of space
  * 
  */
-template <typename IndexVector, typename ScalarVector, typename SegRepType, typename RepfnzType, typename DenseType>
-int LU_copy_to_ucol(const int jcol, const int nseg, SegRepType& segrep, RepfnzType& repfnz ,IndexVector& perm_r, DenseType& dense, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
-{ 
-  typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar; 
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, BlockIndexVector& repfnz ,IndexVector& perm_r, BlockScalarVector& dense, GlobalLU_t& glu)
+{  
   Index ksub, krep, ksupno; 
     
   Index jsupno = glu.supno(jcol);
diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
index 1bda70aaf..6d3271aff 100644
--- a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
@@ -38,8 +38,8 @@
  * \param descendants Number of descendants of each node in the etree
  * \param relax_end last column in a supernode
  */
-template <typename IndexVector>
-void LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_heap_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
 {
   
   // The etree may not be postordered, but its heap ordered  
diff --git a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h
index d5cad49b1..b15ff9c50 100644
--- a/Eigen/src/SparseLU/SparseLU_kernel_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_kernel_bmod.h
@@ -30,7 +30,7 @@ template <int SegSizeAtCompileTime> struct LU_kernel_bmod
   template <typename BlockScalarVector, typename ScalarVector, typename IndexVector>
   EIGEN_DONT_INLINE static void run(const int segsize, BlockScalarVector& dense, ScalarVector& tempv, ScalarVector& lusup, int& luptr, const int nsupr, const int nrow, IndexVector& lsub, const int lptr, const int no_zeros)
   {
-    typedef typename ScalarVector::Scalar Scalar; 
+    typedef typename ScalarVector::Scalar Scalar;
     // First, copy U[*,j] segment from dense(*) to tempv(*)
     // The result of triangular solve is in tempv[*]; 
       // The result of matric-vector update is in dense[*]
diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
index 1b31cc31a..ceb6c5938 100644
--- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
@@ -34,7 +34,7 @@
  * \brief Performs numeric block updates (sup-panel) in topological order.
  * 
  * Before entering this routine, the original nonzeros in the panel
- * were already copied i nto the spa[m,w] ... FIXME to be checked
+ * were already copied i nto the spa[m,w]
  * 
  * \param m number of rows in the matrix
  * \param w Panel size
@@ -48,10 +48,9 @@
  * 
  * 
  */
-template <typename DenseIndexBlock, typename IndexVector, typename ScalarVector>
-void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, DenseIndexBlock& segrep, DenseIndexBlock& repfnz, LU_perfvalues& perfv, LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, LU_perfvalues& perfv, GlobalLU_t& glu)
 {
-  typedef typename ScalarVector::Scalar Scalar; 
   
   int ksub,jj,nextl_col; 
   int fsupc, nsupc, nsupr, nrow; 
@@ -190,9 +189,6 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca
         segsize = krep - kfnz + 1;
         luptr = glu.xlusup(fsupc);    
         
-        // NOTE : Unlike the original implementation in SuperLU, 
-        // there is no update feature for  col-col, 2col-col ... 
-        
         // Perform a trianglar solve and block update, 
         // then scatter the result of sup-col update to dense[]
         no_zeros = kfnz - fsupc; 
@@ -205,4 +201,4 @@ void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, Sca
     
   } // End for each updating supernode
 }
-#endif
\ No newline at end of file
+#endif
diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
index 3581f6d9c..164417897 100644
--- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
@@ -29,12 +29,12 @@
  */
 #ifndef SPARSELU_PANEL_DFS_H
 #define SPARSELU_PANEL_DFS_H
-
-template <typename ScalarVector, typename IndexVector, typename MarkerType, typename Traits>
-void LU_dfs_kernel(const int jj, IndexVector& perm_r,
+template <typename Scalar, typename Index>
+template <typename RepfnzType, typename MarkerType,typename Traits>
+void SparseLUBase<Scalar,Index>::LU_dfs_kernel(const int jj, IndexVector& perm_r,
                    int& nseg, IndexVector& panel_lsub, IndexVector& segrep,
-                   VectorBlock<IndexVector>& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent,
-                   IndexVector& xplore, LU_GlobalLU_t<IndexVector, ScalarVector>& glu,
+                   RepfnzType& repfnz_col, IndexVector& xprune, MarkerType& marker, IndexVector& parent,
+                   IndexVector& xplore, GlobalLU_t& glu,
                    int& nextl_col, int krow, Traits& traits
                   )
 {
@@ -207,8 +207,8 @@ struct LU_panel_dfs_traits
   Index* m_marker;
 };
 
-template <typename MatrixType, typename ScalarVector, typename IndexVector>
-void LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, IndexVector& perm_r, int& nseg, ScalarVector& dense, IndexVector& panel_lsub, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, GlobalLU_t& glu)
 {
   int nextl_col; // Next available position in panel_lsub[*,jj] 
   
diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h
index 4ad49adee..c4a9f1c74 100644
--- a/Eigen/src/SparseLU/SparseLU_pivotL.h
+++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
@@ -52,12 +52,9 @@
  * \return 0 if success, i > 0 if U(i,i) is exactly zero 
  * 
  */
-template <typename IndexVector, typename ScalarVector>
-int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, GlobalLU_t& glu)
 {
-  typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar;
-  typedef typename ScalarVector::RealScalar RealScalar;  
   
   Index fsupc = (glu.xsup)((glu.supno)(jcol)); // First column in the supernode containing the column jcol
   Index nsupc = jcol - fsupc; // Number of columns in the supernode portion, excluding jcol; nsupc >=0
diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h
index f29285bd4..d8c58e039 100644
--- a/Eigen/src/SparseLU/SparseLU_pruneL.h
+++ b/Eigen/src/SparseLU/SparseLU_pruneL.h
@@ -46,12 +46,9 @@
  * \param glu Global LU data
  * 
  */
-template <typename IndexVector, typename ScalarVector, typename BlockIndexVector>
-void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, BlockIndexVector& repfnz, IndexVector& xprune, GlobalLU_t& glu)
 {
-  typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar; 
-    
   // For each supernode-rep irep in U(*,j]
   int jsupno = glu.supno(jcol); 
   int i,irep,irep1; 
diff --git a/Eigen/src/SparseLU/SparseLU_relax_snode.h b/Eigen/src/SparseLU/SparseLU_relax_snode.h
index a9a0a00c1..8db8619c1 100644
--- a/Eigen/src/SparseLU/SparseLU_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_relax_snode.h
@@ -37,8 +37,8 @@
  * \param descendants Number of descendants of each node in the etree
  * \param relax_end last column in a supernode
  */
-template <typename IndexVector>
-void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
+template <typename Scalar, typename Index>
+void SparseLUBase<Scalar,Index>::LU_relax_snode (const int n, IndexVector& et, const int relax_columns, IndexVector& descendants, IndexVector& relax_end)
 {
   
   // compute the number of descendants of each node in the etree
diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
index 18e6a93d2..beea71e31 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -29,11 +29,9 @@
  */
 #ifndef SPARSELU_SNODE_BMOD_H
 #define SPARSELU_SNODE_BMOD_H
-template <typename IndexVector, typename ScalarVector>
-int LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, LU_GlobalLU_t<IndexVector,ScalarVector>& glu)
-{
-  typedef typename ScalarVector::Scalar Scalar; 
-  
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_snode_bmod (const int jcol, const int fsupc, ScalarVector& dense, GlobalLU_t& glu)
+{  
   /* lsub : Compressed row subscripts of ( rectangular supernodes )
    * xlsub :  xlsub[j] is the starting location of the j-th column in lsub(*)
    * lusup : Numerical values of the rectangular supernodes
diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
index edb927cdc..199436cd7 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
@@ -42,55 +42,54 @@
  * \param marker (in/out) working vector
  * \return 0 on success, > 0 size of the memory when memory allocation failed
  */
-  template <typename MatrixType, typename IndexVector, typename ScalarVector>
-  int LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat,  IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+template <typename Scalar, typename Index>
+int SparseLUBase<Scalar,Index>::LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat,  IndexVector& xprune, IndexVector& marker, GlobalLU_t& glu)
+{
+  int mem; 
+  Index nsuper = ++glu.supno(jcol); // Next available supernode number
+  int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
+  int krow,kmark; 
+  for (int i = jcol; i <=kcol; i++)
   {
-    typedef typename IndexVector::Scalar Index; 
-    int mem; 
-    Index nsuper = ++glu.supno(jcol); // Next available supernode number
-    int nextl = glu.xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
-    int krow,kmark; 
-    for (int i = jcol; i <=kcol; i++)
+    // For each nonzero in A(*,i)
+    for (typename MatrixType::InnerIterator it(mat, i); it; ++it)
     {
-      // For each nonzero in A(*,i)
-      for (typename MatrixType::InnerIterator it(mat, i); it; ++it)
+      krow = it.row(); 
+      kmark = marker(krow);
+      if ( kmark != kcol )
       {
-        krow = it.row(); 
-        kmark = marker(krow);
-        if ( kmark != kcol )
+        // First time to visit krow
+        marker(krow) = kcol; 
+        glu.lsub(nextl++) = krow; 
+        if( nextl >= glu.nzlmax )
         {
-          // First time to visit krow
-          marker(krow) = kcol; 
-          glu.lsub(nextl++) = krow; 
-          if( nextl >= glu.nzlmax )
-          {
-            mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
-            if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
-          }
+          mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
+          if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
         }
       }
-      glu.supno(i) = nsuper;
     }
-    
-    // If supernode > 1, then make a copy of the subscripts for pruning
-    if (jcol < kcol)
+    glu.supno(i) = nsuper;
+  }
+  
+  // If supernode > 1, then make a copy of the subscripts for pruning
+  if (jcol < kcol)
+  {
+    Index new_next = nextl + (nextl - glu.xlsub(jcol));
+    while (new_next > glu.nzlmax)
     {
-      Index new_next = nextl + (nextl - glu.xlsub(jcol));
-      while (new_next > glu.nzlmax)
-      {
-        mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
-        if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
-      }
-      Index ifrom, ito = nextl; 
-      for (ifrom = glu.xlsub(jcol); ifrom < nextl;)
-        glu.lsub(ito++) = glu.lsub(ifrom++);
-      for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl;
-      nextl = ito;
+      mem = LUMemXpand<IndexVector>(glu.lsub, glu.nzlmax, nextl, LSUB, glu.num_expansions);
+      if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
     }
-    glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode
-    glu.supno(kcol+1) = nsuper;
-    xprune(kcol) = nextl;
-    glu.xlsub(kcol+1) = nextl;
-    return 0;
+    Index ifrom, ito = nextl; 
+    for (ifrom = glu.xlsub(jcol); ifrom < nextl;)
+      glu.lsub(ito++) = glu.lsub(ifrom++);
+    for (int i = jcol+1; i <=kcol; i++) glu.xlsub(i) = nextl;
+    nextl = ito;
   }
+  glu.xsup(nsuper+1) = kcol + 1; // Start of next available supernode
+  glu.supno(kcol+1) = nsuper;
+  xprune(kcol) = nextl;
+  glu.xlsub(kcol+1) = nextl;
+  return 0;
+}
 #endif
\ No newline at end of file