Checking Data structures and function prototypes

10 files changed, 219 insertions, 189 deletions

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index a4b4fa98b..36b1ce570 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -27,19 +27,12 @@
 #define EIGEN_SPARSE_LU
 
 namespace Eigen {
-  
-template <typename _MatrixType>
-class SparseLU;
 
-#include <Ordering.h>
+
+// Data structure needed by all routines 
 #include <SparseLU_Structs.h>
-#include <SparseLU_Memory.h>
-#include <SparseLU_Utils.h>
-#include <SuperNodalMatrix.h>
+#include <SparseLU_Matrix.h>
 
-#include <SparseLU_Coletree.h>
-#include <SparseLU_heap_relax_snode.h>
-#include <SparseLU_relax_snode.h>
 /**
  * \ingroup SparseLU_Module
  * \brief Sparse supernodal LU factorization for general matrices
@@ -62,7 +55,7 @@ class SparseLU
     typedef Matrix<Index,Dynamic,1> IndexVector;
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
   public:
-    SparseLU():m_isInitialized(true),m_symmetricmode(false),m_fact(DOFACT),m_diagpivotthresh(1.0)
+    SparseLU():m_isInitialized(true),m_symmetricmode(false),m_diagpivotthresh(1.0)
     {
       initperfvalues(); 
     }
@@ -106,7 +99,7 @@ class SparseLU
     }
       
   
-    /** \returns the solution x of \f$ A x = b \f$ using the current decomposition of A.
+    /** \returns the solution X of \f$ A X = b \f$ using the current decomposition of A.
       *
       * \sa compute()
       */
@@ -122,20 +115,34 @@ class SparseLU
   protected:
     // Functions 
     void initperfvalues(); 
-    template <typename IndexVector, typename ScalarVector>
-    int LU_snode_dfs(const int jcol, const int kcol, const IndexVector* asub, 
-                               const IndexVector* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu);
-    
-  template <typename Index, typename ScalarVector>
+  int LU_snode_dfs(const int jcol, const int kcol, const IndexVector* asub,
+                     const IndexVector* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu);
   int LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index fsupc, 
-                      ScalarVector& dense, ScalarVector& tempv, LU_GlobalLu_t& Glu);
+                      ScalarVector& dense, LU_GlobalLU_t& Glu);
+  int LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, 
+                IndexVector& iperm_c, int& pivrow, GlobalLU_t& Glu); 
+  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& Glu);
+  void LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, 
+                     ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, 
+                     IndexVector& repfnz, LU_GlobalLU_t& glu); 
+  int LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg,  
+                    IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, 
+                    IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu); 
+  int LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, 
+                     IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLU_t& Glu); 
+  int LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz,   
+                      IndexVector& perm_r, ScalarVector& dense, LU_GlobalLU_t& glu);  
+  void LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, 
+                 const IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, GlobalLU_t& Glu)
     
     // Variables 
     mutable ComputationInfo m_info;
     bool m_isInitialized;
     bool m_factorizationIsOk;
     bool m_analysisIsOk;
-    fact_t m_fact; 
     NCMatrix m_mat; // The input (permuted ) matrix 
     SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
     NCMatrix m_Ustore; // The upper triangular matrix
@@ -146,7 +153,8 @@ class SparseLU
     ScalarVector m_work; // Scalar work vector 
     IndexVector m_iwork; //Index work vector 
     static LU_GlobalLU_t m_glu; // persistent data to facilitate multiple factors 
-                               // should be defined as a class member
+                               // FIXME All fields of this struct can be defined separately as class members
+                               
     // SuperLU/SparseLU options 
     bool m_symmetricmode;
     
@@ -179,7 +187,10 @@ void SparseLU::initperfvalues()
   m_fillfactor = 20; 
 }
 
-
+// Functions needed by the anaysis phase
+#include <SparseLU_Coletree.h>
+// Ordering interface
+#include <Ordering.h>
 /** 
  * Compute the column permutation to minimize the fill-in (file amd.c )
  * 
@@ -206,7 +217,7 @@ void SparseLU::analyzePattern(const MatrixType& mat)
   
  
   // Apply the permutation to the column of the input  matrix
-  m_mat = mat * m_perm_c;  //FIXME Check if this is valid, check as well how to permute only the index
+  m_mat = mat * m_perm_c; 
     
   // Compute the column elimination tree of the permuted matrix 
   if (m_etree.size() == 0)  m_etree.resize(m_mat.cols());
@@ -234,6 +245,21 @@ void SparseLU::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_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>
+
 /** 
  *  - Numerical factorization 
  *  - Interleaved with the symbolic factorization 
@@ -284,7 +310,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
   idx += m;
   VectorBlock<IndexVector> xplore(m_iwork, idx, m); 
   idx += m; 
-  VectorBlock<IndexVector> repnfnz(m_iwork, idx, maxpanel);
+  VectorBlock<IndexVector> repfnz(m_iwork, idx, maxpanel);
   idx += maxpanel;
   VectorBlock<IndexVector> panel_lsub(m_iwork, idx, maxpanel)
   idx += maxpanel; 
@@ -324,8 +350,6 @@ void  SparseLU::factorize(const MatrixType& matrix)
     
   supno(0) = IND_EMPTY; 
   xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = 0;
-  int panel_size = m_panel_size; 
-  int wdef = m_panel_size; // upper bound on panel width
   
   // 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
@@ -348,9 +372,9 @@ void  SparseLU::factorize(const MatrixType& matrix)
       info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker); 
       if ( info ) 
       {
+        std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n";
         m_info = NumericalIssue; 
         m_factorizationIsOk = false; 
-        std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n";
         return; 
       }
       nextu = xusub(jcol); //starting location of column jcol in ucol
@@ -377,13 +401,14 @@ void  SparseLU::factorize(const MatrixType& matrix)
           dense(it.row()) = it.val();
         
         // Numeric update within the snode 
-        LU_snode_bmod(icol, jsupno, fsupc, dense, tempv); 
+        LU_snode_bmod(icol, jsupno, fsupc, dense, glu); 
         
         // Eliminate the current column 
         info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r, m_iperm_c, pivrow, m_glu); 
         if ( info ) 
         {
           m_info = NumericalIssue; 
+          std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl; 
           m_factorizationIsOk = false; 
           return; 
         }
@@ -394,7 +419,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
     { // Work on one panel of panel_size columns
       
       // Adjust panel size so that a panel won't overlap with the next relaxed snode. 
-      panel_size = w_def;
+      int panel_size = wdef; // upper bound on panel width
       for (k = jcol + 1; k < std::min(jcol+panel_size, n); k++)
       {
         if (relax_end(k) != IND_EMPTY) 
@@ -419,31 +444,33 @@ void  SparseLU::factorize(const MatrixType& matrix)
         
         nseg = nseg1; // begin after all the panel segments
         //Depth-first-search for the current column
-        VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); //FIXME
-        VectorBlock<IndexVector> repfnz_k(repfnz, k, m); //FIXME 
+        VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m);
+        VectorBlock<IndexVector> repfnz_k(repfnz, k, m);
         info = LU_column_dfs(m, jj, perm_r, nseg, panel_lsub(k), segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 
         if ( !info ) 
         {
+          std::cerr << "UNABLE TO EXPAND MEMORY IN COLUMN_DFS() \n";
           m_info = NumericalIssue; 
           m_factorizationIsOk = false; 
           return; 
         }
         // Numeric updates to this column 
-        VectorBlock<IndexVector> dense_k(dense, k, m); //FIXME 
-        VectorBlock<IndexVector> segrep_k(segrep, nseg1, m) // FIXME Check the length
+        VectorBlock<IndexVector> dense_k(dense, k, m); 
+        VectorBlock<IndexVector> segrep_k(segrep, nseg1, m); 
         info = 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";
           m_info = NumericalIssue; 
           m_factorizationIsOk = false; 
           return; 
         }
         
         // Copy the U-segments to ucol(*)
-        //FIXME Check that repfnz_k, dense_k... have stored references to modified columns
         info = LU_copy_to_col(jj, nseg, segrep, repfnz_k, perm_r, dense_k, m_glu); 
         if ( info ) 
         {
+          std::cerr << "UNABLE TO EXPAND MEMORY IN COPY_TO_UCOL() \n";
           m_info = NumericalIssue; 
           m_factorizationIsOk = false; 
           return; 
@@ -453,6 +480,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r, iperm_c, pivrow, m_glu);
         if ( info ) 
         {
+          std::cerr<< "THE MATRIX IS STRUCTURALLY SINGULAR ... ZERO COLUMN AT " << info <<std::endl; 
           m_info = NumericalIssue; 
           m_factorizationIsOk = false; 
           return; 
@@ -472,7 +500,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
     } // end else 
   } // end for -- end elimination 
   
-  // Adjust row permutation in the case of rectangular matrices
+  // Adjust row permutation in the case of rectangular matrices... Deprecated
   if (m > n ) 
   {
     k = 0; 
@@ -504,18 +532,18 @@ void  SparseLU::factorize(const MatrixType& matrix)
 }
 
 template<typename Rhs, typename Dest>
-bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
+bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &X) const
 {
   eigen_assert(m_isInitialized && "The matrix should be factorized first");
   EIGEN_STATIC_ASSERT((Dest::Flags&RowMajorBit)==0,
                      THIS_METHOD_IS_ONLY_FOR_COLUMN_MAJOR_MATRICES);
   
-  x = b; /* on return, x is overwritten by the computed solution */
+  X = b; /* on return, X is overwritten by the computed solution */
   
   int nrhs = b.cols(); 
   
   // Permute the right hand side to form Pr*B
-  x = m_perm_r * x; 
+  X = m_perm_r * X; 
   
   // Forward solve PLy = Pb; 
   Index fsupc; // First column of the current supernode 
@@ -547,7 +575,7 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
         {
           irow = m_Lstore.rowIndex()[iptr]; 
           ++luptr; 
-          x(irow, j) -= x(fsupc, j) * Lval[luptr]; 
+          X(irow, j) -= X(fsupc, j) * Lval[luptr]; 
         }
       }
     }
@@ -558,8 +586,8 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
       // Triangular solve 
       luptr = m_Lstore.colIndexPtr()[fsupc]; //FIXME Should be outside the loop 
       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); 
-//       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride > u( &(x(fsupc,0)), nsupc, nrhs, OuterStride<>(x.rows()) );
-      Matrix<Scalar,Dynamic,Dynamic>& u = x.block(fsupc, 0, nsupc, nrhs); //FIXME Check this 
+//       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride > u( &(X(fsupc,0)), nsupc, nrhs, OuterStride<>(X.rows()) );
+      Matrix<Scalar,Dynamic,Dynamic>& u = X.block(fsupc, 0, nsupc, nrhs); //FIXME Check this 
       u = A.triangularView<Lower>().solve(u); 
       
       // Matrix-vector product 
@@ -573,7 +601,7 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
         for (i = 0; i < nrow; i++)
         {
           irow = m_Lstore.rowIndex()[iptr]; 
-          x(irow, j) -= work(i, j); // Scatter operation
+          X(irow, j) -= work(i, j); // Scatter operation
           work(i, j) = Scalar(0); 
           iptr++;
         }
@@ -594,13 +622,13 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
     {
       for (j = 0; j < nrhs; j++)
       {
-	x(fsupc, j) /= Lval[luptr]; 
+        X(fsupc, j) /= Lval[luptr]; 
       }
     }
     else 
     {
       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(Lval[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); 
-      Matrix<Scalar,Dynamic,Dynamic>& u = x.block(fsupc, 0, nsupc, nrhs); 
+      Matrix<Scalar,Dynamic,Dynamic>& u = X.block(fsupc, 0, nsupc, nrhs); 
       u = A.triangularView<Upper>().solve(u); 
     }
     
@@ -608,17 +636,17 @@ bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &x) const
     {
       for (jcol = fsupc; jcol < fsupc + nsupc; jcol++)
       {
-	for (i = m_Ustore.outerIndexPtr()[jcol]; i < m_Ustore.outerIndexPtr()[jcol]; i++)
-	{
-	  irow = m_Ustore.InnerIndices()[i]; 
-	  x(irow, j) -= x(irow, jcol) * m_Ustore.Values()[i];
-	}
+        for (i = m_Ustore.outerIndexPtr()[jcol]; i < m_Ustore.outerIndexPtr()[jcol]; i++)
+        {
+          irow = m_Ustore.InnerIndices()[i]; 
+          X(irow, j) -= X(irow, jcol) * m_Ustore.Values()[i];
+        }
       }
     }
   } // End For U-solve
   
   // Permute back the solution 
-  x = x * m_perm_c; 
+  X = m_perm_c * X; 
   
   return true; 
 }
diff --git a/Eigen/src/SparseLU/SparseLU_Structs.h b/Eigen/src/SparseLU/SparseLU_Structs.h
index 1394eccdf..618d05eac 100644
--- a/Eigen/src/SparseLU/SparseLU_Structs.h
+++ b/Eigen/src/SparseLU/SparseLU_Structs.h
@@ -82,19 +82,12 @@
  */
 #ifndef EIGEN_LU_STRUCTS
 #define EIGEN_LU_STRUCTS
-namespace Eigen {
-
-#define LU_NBR_MEMTYPE 4 /* 0: lusup
-                        1: ucol 
-                        2: lsub 
-                        3: usub */
-typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD, MY_PERMC} colperm_t; 
-typedef enum {DOFACT, SamePattern, Factored} fact_t; 
 typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} MemType; 
 
 
 template <typename ScalarVector, typename IndexVector>
 struct {
+  typedef typename IndexVector::Index Index; 
   IndexVector xsup; //First supernode column ... xsup(s) points to the beginning of the s-th supernode
   IndexVector supno; // Supernode number corresponding to this column (column to supernode mapping)
   ScalarVector  lusup; // nonzero values of L ordered by columns 
@@ -113,5 +106,4 @@ struct {
   int   num_expansions; 
 } GlobalLU_t;
 
-}// End namespace Eigen 
 #endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_column_bmod.h b/Eigen/src/SparseLU/SparseLU_column_bmod.h
index bed4f9519..965a0c0ad 100644
--- a/Eigen/src/SparseLU/SparseLU_column_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_column_bmod.h
@@ -54,35 +54,40 @@
  * \param segrep segment representative ...
  * \param repfnz ??? First nonzero column in each row ???  ...
  * \param fpanelc First column in the current panel
- * \param Glu Global LU data. 
+ * \param glu Global LU data. 
  * \return 0 - successful return 
  *         > 0 - number of bytes allocated when run out of space
  * 
  */
 template <typename ScalarVector, typename IndexVector>
-int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLu_t& Glu)
+int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, int fpanelc, LU_GlobalLU_t& glu)
 {
     
-  int jsupno, k, ksub, krep, krep_ind, ksupno; 
+  int  jsupno, k, ksub, krep, krep_ind, ksupno; 
+  int fsupc, nsupc, nsupr, luptr, kfnz, no_zeros; 
   /* krep = representative of current k-th supernode
     * fsupc =  first supernodal column
     * nsupc = number of columns in a supernode
     * nsupr = number of rows in a supernode
     * luptr = location of supernodal LU-block in storage
     * kfnz = first nonz in the k-th supernodal segment
-    * no-zeros = no lf leading zeros in a supernodal U-segment
+    * no_zeros = no lf leading zeros in a supernodal U-segment
     */
-  IndexVector& xsup = Glu.xsup; 
-  IndexVector& supno = Glu.supno; 
-  IndexVector& lsub = Glu.lsub; 
-  IndexVector& xlsub = Glu.xlsub; 
-  IndexVector& xlusup = Glu.xlusup; 
-  ScalarVector& lusup = Glu.lusup; 
-  Index& nzlumax = Glu.nzlumax; 
+  IndexVector& xsup = glu.xsup; 
+  IndexVector& supno = glu.supno; 
+  IndexVector& lsub = glu.lsub; 
+  IndexVector& xlsub = glu.xlsub; 
+  IndexVector& xlusup = glu.xlusup; 
+  ScalarVector& lusup = glu.lusup; 
+  Index& nzlumax = glu.nzlumax; 
   
   int jsupno = supno(jcol);
   // For each nonzero supernode segment of U[*,j] in topological order 
   k = nseg - 1; 
+  int d_fsupc; // distance between the first column of the current panel and the 
+               // first column of the current snode
+  int fst_col; // First column within small LU update
+  int segsize; 
   for (ksub = 0; ksub < nseg; ksub++)
   {
     krep = segrep(k); k--; 
@@ -110,35 +115,36 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense
       krep_ind = lptr + nsupc - 1; 
       
       // NOTE  Unlike the original implementation in SuperLU, the only feature  
-      //  here is a sup-col update. 
+      // 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; 
       // First, copy U[*,j] segment from dense(*) to tempv(*)
       isub = lptr + no_zeros; 
-      for (i = 0; i ww segsize; i++)
+      for (i = 0; i < segsize; i++)
       {
         irow = lsub(isub); 
-        tempv(i) = densee(irow); 
+        tempv(i) = dense(irow); 
         ++isub; 
       }
       // Dense triangular solve -- start effective triangle
       luptr += nsupr * no_zeros + no_zeros; 
       // Form Eigen matrix and vector 
       Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) );
-      Map<ScalarVector> u(tempv.data(), segsize);
+      VectorBlock<ScalarVector> u(tempv, 0, segsize);
+      
       u = A.triangularView<Lower>().solve(u); 
       
       // Dense matrix-vector product y <-- A*x 
       luptr += segsize; 
-      new (&A) (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); 
-      Map<ScalarVector> l( &(tempv.data()[segsize]), segsize); 
+      new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); 
+      VectorBlock<ScalarVector> l(tempv, segsize, nrow); 
       l= A * u;
       
       // Scatter tempv[] into SPA dense[] as a temporary storage 
       isub = lptr + no_zeros; 
-      for (i = 0; i w segsize; i++)
+      for (i = 0; i < segsize; i++)
       {
         irow = lsub(isub); 
         dense(irow) = tempv(i); 
@@ -150,8 +156,8 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense
       for (i = 0; i < nrow; i++)
       {
         irow = lsub(isub); 
-        dense(irow) -= tempv(segsize + i); 
-        tempv(segsize + i) = Scalar(0.0); 
+        dense(irow) -= l(i); 
+        l(i) = Scalar(0.0); 
         ++isub; 
       }
     } // end if jsupno 
@@ -165,9 +171,9 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense
   new_next = nextlu + xlsub(fsupc + 1) - xlsub(fsupc); 
   while (new_next > nzlumax )
   {
-    mem = LUmemXpand<Scalar>(Glu.lusup, nzlumax, nextlu, LUSUP, Glu);  
+    mem = LUmemXpand<Scalar>(glu.lusup, nzlumax, nextlu, LUSUP, glu);  
     if (mem) return mem; 
-    lsub = Glu.lsub; //FIXME Why is it updated here. 
+    //lsub = glu.lsub; // Should not be updated here 
   }
   
   for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
@@ -183,8 +189,8 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense
   /* For more updates within the panel (also within the current supernode),
    * should start from the first column of the panel, or the first column
    * of the supernode, whichever is bigger. There are two cases:
-   *  1) fsupc < fpanelc, then fst_col <- fpanelc
-   *  2) fsupc >= fpanelc, then fst_col <-fsupc
+   *  1) fsupc < fpanelc, then fst_col <-- fpanelc
+   *  2) fsupc >= fpanelc, then fst_col <-- fsupc
    */
   fst_col = std::max(fsupc, fpanelc); 
   
@@ -203,11 +209,11 @@ int SparseLU::LU_column_bmod(const int jcol, const int nseg, ScalarVector& dense
     // points to the beginning of jcol in snode L\U(jsupno) 
     ufirst = xlusup(jcol) + d_fsupc; 
     Map<Matrix<Scalar,Dynamic,Dynamic>, 0,  OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); 
-    Map<ScalarVector> l( &(lusup.data()[ufirst]), nsupc ); 
+    VectorBlock<ScalarVector> u(lusup, ufirst, nsupc); 
     u = A.triangularView().solve(u); 
     
     new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); 
-    Map<ScalarVector> l( &(lusup.data()[ufirst+nsupc]), nsupr ); 
+    VectorBlock<ScalarVector> l(lusup, ufirst+nsupc, nrow); 
     l = l - A * u;
     
   } // End if fst_col
diff --git a/Eigen/src/SparseLU/SparseLU_column_dfs.h b/Eigen/src/SparseLU/SparseLU_column_dfs.h
index 1c832d60e..7fda536a9 100644
--- a/Eigen/src/SparseLU/SparseLU_column_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_column_dfs.h
@@ -65,13 +65,13 @@
  * \param marker
  * \param parent
  * \param xplore
- * \param Glu global LU data 
+ * \param glu global LU data 
  * \return 0 success
  *         > 0 number of bytes allocated when run out of space
  * 
  */
-template <typename IndexVector>
-int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg  IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLu_t& Glu)
+template <typename IndexVector, typename ScalarVector>
+int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, IndexVector& nseg  IndexVector& lsub_col, IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, IndexVector& marker, IndexVector& parent, IndexVector& xplore, LU_GlobalLU_t& glu)
 {
   typedef typename IndexVector::IndexVector; 
   
@@ -85,15 +85,15 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In
   int xdfs, maxdfs, kpar;
   int mem; 
   // Initialize pointers 
-  IndexVector& xsup = Glu.xsup; 
-  IndexVector& supno = Glu.supno; 
-  IndexVector& lsub = Glu.lsub; 
-  IndexVector& xlsub = Glu.xlsub; 
-  IndexVector& nzlmax = Glu.nzlmax; 
+  IndexVector& xsup = glu.xsup; 
+  IndexVector& supno = glu.supno; 
+  IndexVector& lsub = glu.lsub; 
+  IndexVector& xlsub = glu.xlsub; 
+  IndexVector& nzlmax = glu.nzlmax; 
   
   nsuper = supno(jcol); 
   jsuper = nsuper; 
-  nextl = xlsup(jcol); 
+  nextl = xlsub(jcol); 
   VectorBlock<IndexVector> marker2(marker, 2*m, m); 
   // For each nonzero in A(*,jcol) do dfs 
   for (k = 0; lsub_col[k] != IND_EMPTY; k++) 
@@ -106,16 +106,16 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In
     if (kmark == jcol) continue; 
     
     // For each unmarker nbr krow of jcol
-    //  krow is in L: place it in structure of L(*,jcol)
     marker2(krow) = jcol; 
     kperm = perm_r(krow); 
     
     if (kperm == IND_EMPTY ) 
     {
+      //  krow is in L: place it in structure of L(*,jcol)
       lsub(nextl++) = krow; // krow is indexed into A
       if ( nextl >= nzlmax )
       {
-        mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, Glu); 
+        mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu); 
         if ( mem ) return mem; 
       }
       if (kmark != jcolm1) jsuper = IND_EMPTY; // Row index subset testing
@@ -157,13 +157,13 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In
               marker2(kchild) = jcol; 
               chperm = perm_r(kchild); 
               
-              // if kchild is in L: place it in L(*,k)
               if (chperm == IND_EMPTY)
               {
+                // if kchild is in L: place it in L(*,k)
                 lsub(nextl++) = kchild; 
                 if (nextl >= nzlmax)
                 {
-                   mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB); 
+                   mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu); 
                    if (mem) return mem; 
                 }
                 if (chmark != jcolm1) jsuper = IND_EMPTY; 
@@ -201,7 +201,7 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In
           // place supernode-rep krep in postorder DFS.
           // backtrack dfs to its parent
           
-          segrep(nseg) = ;krep; 
+          segrep(nseg) = krep; 
           ++nseg; 
           kpar = parent(krep); // Pop from stack, mimic recursion
           if (kpar == IND_EMPTY) break; // dfs done 
@@ -217,7 +217,7 @@ int SparseLU::LU_column_dfs(const int m, const int jcol, IndexVector& perm_r, In
     
   } // for each nonzero ... 
   
-  // check to see if j belongs in the same supeprnode as j-1
+  // check to see if j belongs in the same supernode as j-1
   if ( jcol == 0 )
   { // Do nothing for column 0 
     nsuper = supno(0) = 0 ;
diff --git a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
index dc53edcfb..c97bc6aa4 100644
--- a/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
+++ b/Eigen/src/SparseLU/SparseLU_copy_to_ucol.h
@@ -50,28 +50,28 @@
  * \param jcol current column to update
  * \param nseg Number of segments in the U part
  * \param segrep segment representative ...
- * \param repfnz ??? First nonzero column in each row ???  ...
+ * \param repfnz First nonzero column in each row  ...
  * \param perm_r Row permutation 
  * \param dense Store the full representation of the column
- * \param Glu Global LU data. 
+ * \param glu Global LU data. 
  * \return 0 - successful return 
  *         > 0 - number of bytes allocated when run out of space
  * 
  */
 template <typename ScalarVector, typename IndexVector>
-int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz, IndexVector& perm_r, ScalarVector& dense, LU_GlobalLu_t& Glu)
+int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segrep, IndexVector& repfnz, IndexVector& perm_r, ScalarVector& dense, LU_GlobalLU_t& glu)
 { 
   Index ksupno, k, ksub, krep, ksupno; 
   typedef typename IndexVector::Index; 
   
-  IndexVector& xsup = Glu.xsup; 
-  IndexVector& supno = Glu.supno; 
-  IndexVector& lsub = Glu.lsub; 
-  IndexVector& xlsub = Glu.xlsub; 
+  IndexVector& xsup = glu.xsup; 
+  IndexVector& supno = glu.supno; 
+  IndexVector& lsub = glu.lsub; 
+  IndexVector& xlsub = glu.xlsub; 
   ScalarVector& ucol = GLu.ucol; 
-  IndexVector& usub = Glu.usub; 
-  IndexVector& xusub = Glu.xusub;
-  Index& nzumax = Glu.nzumax; 
+  IndexVector& usub = glu.usub; 
+  IndexVector& xusub = glu.xusub;
+  Index& nzumax = glu.nzumax; 
   
   Index jsupno = supno(jcol);
   
@@ -95,12 +95,11 @@ int SparseLU::LU_copy_to_ucol(const int jcol, const int nseg, IndexVector& segre
         new_next = nextu + segsize; 
         while (new_next > nzumax) 
         {
-          mem = LU_MemXpand<ScalarVector>(ucol, nzumax, nextu, UCOL, Glu); 
+          mem = LU_MemXpand<ScalarVector>(ucol, nzumax, nextu, UCOL, glu); 
           if (mem) return mem; 
-          mem = LU_MemXpand<Index>(usub, nzumax, nextu, USUB, Glu); 
+          mem = LU_MemXpand<Index>(usub, nzumax, nextu, USUB, glu); 
           if (mem) return mem; 
           
-          lsub = Glu.lsub; //FIXME Why setting this as well ??
         }
         
         for (i = 0; i < segsize; i++)
diff --git a/Eigen/src/SparseLU/SparseLU_panel_bmod.h b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
index 93daa938c..212ecfa6a 100644
--- a/Eigen/src/SparseLU/SparseLU_panel_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_panel_bmod.h
@@ -56,21 +56,21 @@
  * \param nseg Number of segments in the U part
  * \param dense Store the full representation of the panel 
  * \param tempv working array 
- * \param segrep in ...
- * \param repfnz in ...
- * \param Glu Global LU data. 
+ * \param segrep segment representative... first row in the segment
+ * \param repfnz First nonzero rows
+ * \param glu Global LU data. 
  * 
  * 
  */
-template <typename VectorType>
-void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, VectorType& dense, VectorType& tempv, VectorXi& segrep, VectorXi& repfnz, LU_GlobalLu_t& Glu)
+template <typename IndexVector, typename ScalarVector>
+void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int nseg, ScalarVector& dense, ScalarVector& tempv, IndexVector& segrep, IndexVector& repfnz, LU_GlobalLU_t& glu)
 {
-  VectorXi& xsup = Glu.xsup; 
-  VectorXi& supno = Glu.supno; 
-  VectorXi& lsub = Glu.lsub; 
-  VectorXi& xlsub = Glu.xlsub; 
-  VectorXi& xlusup = Glu.xlusup; 
-  VectorType& lusup = Glu.lusup; 
+  IndexVector& xsup = glu.xsup; 
+  IndexVector& supno = glu.supno; 
+  IndexVector& lsub = glu.lsub; 
+  IndexVector& xlsub = glu.xlsub; 
+  IndexVector& xlusup = glu.xlusup; 
+  ScalarVector& lusup = glu.lusup; 
   
   int i,ksub,jj,nextl_col,irow; 
   int fsupc, nsupc, nsupr, nrow; 
@@ -96,10 +96,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
     nrow = nsupr - nsupc; 
     lptr = xlsub(fsupc); 
     krep_ind = lptr + nsupc - 1; 
-    
-    repfnz_col = repfnz; 
-    dense_col = dense; 
-    
+        
     // NOTE : Unlike the original implementation in SuperLU, the present implementation
     // does not include a 2-D block update. 
     
@@ -107,8 +104,8 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
     for (jj = jcol; jj < jcol + w; jj++)
     {
       nextl_col = (jj-jcol) * m; 
-      VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero column index for each row
-      VectorBLock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here
+      VectorBlock<IndexVector> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero column index for each row
+      VectorBLock<IndexVector> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here
       
       kfnz = repfnz_col(krep); 
       if ( kfnz == IND_EMPTY ) 
@@ -123,8 +120,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
       // Perform a trianglar solve and block update, 
       // then scatter the result of sup-col update to dense[]
       no_zeros = kfnz - fsupc; 
-      
-      // Copy U[*,j] segment from dense[*] to tempv[*] :
+      // 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_col[*]
       isub = lptr + no_zeros; 
@@ -138,19 +134,21 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
       luptr += nsupr * no_zeros + no_zeros; 
       // triangular solve with Eigen
       Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); 
-      Map<Matrix<Scalar,Dynamic,1> > u( tempv.data(), segsize); 
+//       Map<Matrix<Scalar,Dynamic,1> > u( tempv.data(), segsize); 
+      VectorBlock<ScalarVector> u(tempv, 0, segsize);
       u = A.triangularView<Lower>().solve(u); 
       
       luptr += segsize; 
       // Dense Matrix vector product y <-- A*x; 
       new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); 
-      Map<VectorType> l( &(tempv.data()[segsize]), segsize); 
+//       Map<ScalarVector> l( &(tempv.data()[segsize]), nrow); 
+      VectorBlock<ScalarVector> l(tempv, segsize, nrow); 
       l= A * u;
       
       // Scatter tempv(*) into SPA dense(*) such that tempv(*) 
       // can be used for the triangular solve of the next 
       // column of the panel. The y will be copied into ucol(*) 
-      // after the whole panel has been finished.
+      // after the whole panel has been finished... after column_dfs() and column_bmod()
       
       isub = lptr + no_zeros; 
       for (i = 0; i < segsize; i++) 
@@ -166,8 +164,8 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
      for (i = 0; i < nrow; i++) 
      {
        irow = lsub(isub); 
-       dense_col(irow) -= tempv(segsize + i); 
-       tempv(segsize + i) = 0; 
+       dense_col(irow) -= l(i); 
+       l(i) = Scalar(0); 
        ++isub; 
      }
      
diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
index 97e5121db..d3c2906b2 100644
--- a/Eigen/src/SparseLU/SparseLU_panel_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
@@ -62,45 +62,50 @@
  *    marker[i] == jj, if i was visited during dfs of current column jj;
  *    marker1[i] >= jcol, if i was visited by earlier columns in this panel; 
  * 
- * \param m number of rows in the matrix
- * \param w Panel size
- * \param jcol Starting  column of the panel
- * \param A Input matrix in column-major storage
- * \param perm_r Row permutation
- * \param nseg Number of U segments
- * ...
+ * \param [in]m number of rows in the matrix
+ * \param [in]w Panel size
+ * \param [in]jcol Starting  column of the panel
+ * \param [in]A Input matrix in column-major storage
+ * \param [in]perm_r Row permutation
+ * \param [out]nseg Number of U segments
+ * \param [out]dense Accumulate the column vectors of the panel
+ * \param [out]panel_lsub Subscripts of the row in the panel 
+ * \param [out]segrep Segment representative i.e first nonzero row of each segment
+ * \param [out]repfnz First nonzero location in each row
+ * \param [out]xprune 
+ * \param [out]marker 
+ * 
  * 
  */
-template <typename MatrixType, typename VectorType>
-void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, int& nseg, VectorType& dense,  VectorXi& panel_lsub, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu)
+template <typename MatrixType, typename IndexVector, typename ScalarVector>
+void SparseLU::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& Glu)
 {
   
   int jj; // Index through each column in the panel 
   int nextl_col; // Next available position in panel_lsub[*,jj] 
   int krow; // Row index of the current element 
   int kperm; // permuted row index
-  int krep; // Supernode reprentative of the current row
+  int krep; // Supernode representative of the current row
   int kmark; 
   int chperm, chmark, chrep, oldrep, kchild; 
   int myfnz; // First nonzero element in the current column
   int xdfs, maxdfs, kpar;
   
   // Initialize pointers 
-//   VectorXi& marker1 = marker.block(m, m); 
-  VectorBlock<VectorXi> marker1(marker, m, m); 
+//   IndexVector& marker1 = marker.block(m, m); 
+  VectorBlock<IndexVector> marker1(marker, m, m); 
   nseg = 0; 
-  VectorXi& xsup = Glu.xsup; 
-  VectorXi& supno = Glu.supno; 
-  VectorXi& lsub = Glu.lsub; 
-  VectorXi& xlsub = Glu.xlsub; 
+  IndexVector& xsup = Glu.xsup; 
+  IndexVector& supno = Glu.supno; 
+  IndexVector& lsub = Glu.lsub; 
+  IndexVector& xlsub = Glu.xlsub; 
   // For each column in the panel 
   for (jj = jcol; jj < jcol + w; jj++) 
   {
     nextl_col = (jj - jcol) * m; 
     
-    //FIXME 
-    VectorBlock<VectorXi> repfnz_col(repfnz.segment(nextl_col, m)); // First nonzero location in each row
-    VectorBlock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Accumulate a column vector here
+    VectorBlock<IndexVector> repfnz_col(repfnz, nextl_col, m); // First nonzero location in each row
+    VectorBlock<IndexVector> dense_col(dense,nextl_col, m); // Accumulate a column vector here
     
     
     // For each nnz in A[*, jj] do depth first search
diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h
index 3bfe14e7e..32da92481 100644
--- a/Eigen/src/SparseLU/SparseLU_pivotL.h
+++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
@@ -67,28 +67,30 @@
  * \return 0 if success, i > 0 if U(i,i) is exactly zero 
  * 
  */
-template <typename Scalar>
-int SparseLU::LU_pivotL(const int jcol, const Scalar u, VectorXi& perm_r, VectorXi& iperm_c, int& pivrow, GlobalLU_t& Glu)
+template <typename IndexVector, typename ScalarVector>
+int SparseLU::LU_pivotL(const int jcol, const RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, GlobalLU_t& Glu)
 {
+  typedef typename IndexVector::Index Index; 
+  typedef typename ScalarVector::Scalar Scalar; 
   // Initialize pointers 
-  VectorXi& lsub = Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
-  VectorXi& xlsub = Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
-  Scalar* lusup = Glu.lusup.data(); // Numerical values of the rectangular supernodes
-  VectorXi& xlusup = Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
+  IndexVector& lsub = Glu.lsub; // Compressed row subscripts of L rectangular supernodes.
+  IndexVector& xlsub = Glu.xlsub; // pointers to the beginning of each column subscript in lsub
+  ScalarVector& lusup = Glu.lusup; // Numerical values of L ordered by columns 
+  IndexVector& xlusup = Glu.xlusup; //  pointers to the beginning of each colum in lusup
   
   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
   Index lptr = xlsub(fsupc); // pointer to the starting location of the row subscripts for this supernode portion
   Index nsupr = xlsub(fsupc+1) - lptr; // Number of rows in the supernode
-  Scalar* lu_sup_ptr = &(lusup[xlusup(fsupc)]); // Start of the current supernode
-  Scalar* lu_col_ptr = &(lusup[xlusup(jcol)]); // Start of jcol in the supernode
+  Scalar* lu_sup_ptr = &(lusup.data()[xlusup(fsupc)]); // Start of the current supernode
+  Scalar* lu_col_ptr = &(lusup.data()[xlusup(jcol)]); // Start of jcol in the supernode
   Index* lsub_ptr = &(lsub.data()[lptr]); // Start of row indices of the supernode
   
   // Determine the largest abs numerical value for partial pivoting 
   Index diagind = iperm_c(jcol); // diagonal index 
   Scalar pivmax = 0.0; 
   Index pivptr = nsupc; 
-  Index diag = -1; 
+  Index diag = IND_EMPTY; 
   Index old_pivptr = nsupc; 
   Scalar rtemp;
   for (isub = nsupc; isub < nsupr; ++isub) {
@@ -127,7 +129,7 @@ int SparseLU::LU_pivotL(const int jcol, const Scalar u, VectorXi& perm_r, Vector
   // Interchange row subscripts
   if (pivptr != nsupc )
   {
-    std::swap( lsub_ptr(pivptr), lsub_ptr(nsupc) );
+    std::swap( lsub_ptr[pivptr], lsub_ptr[nsupc] );
     // Interchange numerical values as well, for the two rows in the whole snode
     // such that L is indexed the same way as A
     for (icol = 0; icol <= nsupc; icol++)
diff --git a/Eigen/src/SparseLU/SparseLU_pruneL.h b/Eigen/src/SparseLU/SparseLU_pruneL.h
index 687717d52..dd092b778 100644
--- a/Eigen/src/SparseLU/SparseLU_pruneL.h
+++ b/Eigen/src/SparseLU/SparseLU_pruneL.h
@@ -47,35 +47,35 @@
 /**
  * \brief Prunes the L-structure.
  *
- * It prunes the L-structure  of supernodes whose L-structure constains the current pivot row "pivrow"
+ * It prunes the L-structure  of supernodes whose L-structure contains the current pivot row "pivrow"
  * 
  * 
  * \param jcol The current column of L
  * \param [in]perm_r Row permutation
  * \param [out]pivrow  The pivot row
- * \param nseg Number of segments ???
+ * \param nseg Number of segments
  * \param segrep 
  * \param repfnz
  * \param [out]xprune 
  * \param Glu Global LU data
  * 
  */
-template <typename VectorType>
-void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivrow, const int nseg, const VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, GlobalLU_t& Glu)
+template <typename IndexVector, typename ScalarVector>
+void SparseLU::LU_pruneL(const int jcol, const IndexVector& perm_r, const int pivrow, const int nseg, const IndexVector& segrep, IndexVector& repfnz, IndexVector& xprune, GlobalLU_t& Glu)
 {
   // Initialize pointers 
-  VectorXi& xsup = Glu.xsup; 
-  VectorXi& supno = Glu.supno; 
-  VectorXi& lsub = Glu.lsub; 
-  VectorXi& xlsub = Glu.xlsub; 
-  VectorType& lusup = Glu.lusup;
-  VectorXi& xlusup = Glu.xlusup; 
+  IndexVector& xsup = Glu.xsup; 
+  IndexVector& supno = Glu.supno; 
+  IndexVector& lsub = Glu.lsub; 
+  IndexVector& xlsub = Glu.xlsub; 
+  ScalarVector& lusup = Glu.lusup;
+  IndexVector& xlusup = Glu.xlusup; 
   
   // For each supernode-rep irep in U(*,j]
   int jsupno = supno(jcol); 
   int i,irep,irep1; 
   bool movnum, do_prune = false; 
-  int kmin, kmax, ktemp, minloc, maxloc; 
+  Index kmin, kmax, ktemp, minloc, maxloc; 
   for (i = 0; i < nseg; i++)
   {
     irep = segrep(i); 
@@ -125,9 +125,7 @@ void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivro
           {
             // kmin below pivrow (not yet pivoted), and kmax
             // above pivrow: interchange the two suscripts
-            ktemp = lsub(kmin);
-            lsub(kmin) = lsub(kmax);
-            lsub(kmax) = ktemp; 
+            std::swap(lsub(kmin), lsub(kmax)); 
             
             // If the supernode has only one column, then we 
             // only keep one set of subscripts. For any subscript
@@ -144,7 +142,7 @@ void SparseLU::LU_pruneL(const int jcol, const VectorXi& perm_r, const int pivro
           }
         } // end while 
         
-        xprune(irep) = kmin; 
+        xprune(irep) = kmin;  //Pruning 
       } // end if do_prune 
     } // end pruning 
   } // End for each U-segment
diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
index 6130a5622..1d6bed8bb 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -47,13 +47,13 @@ namespace internal {
 #define SPARSELU_SNODE_BMOD_H
 template <typename Index, typename ScalarVector>
 int SparseLU::LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index fsupc, 
-                              ScalarVector& dense, ScalarVector& tempv, LU_GlobalLu_t& Glu)
+                              ScalarVector& dense, LU_GlobalLU_t& glu)
 {
   typedef typename Matrix<Index, Dynamic, Dynamic> IndexVector;
-  IndexVector& lsub = Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
-  IndexVector& xlsub = Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
-  ScalarVector& lusup = Glu.lusup; // Numerical values of the rectangular supernodes
-  IndexVector& xlusup = Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
+  IndexVector& lsub = glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
+  IndexVector& xlsub = glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
+  ScalarVector& lusup = glu.lusup; // Numerical values of the rectangular supernodes
+  IndexVector& xlusup = glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
   
   int nextlu = xlusup(jcol); // Starting location of the next column to add 
   int irow, isub; 
@@ -75,14 +75,16 @@ int SparseLU::LU_dsnode_bmod (const Index jcol, const Index jsupno, const Index
     
     int nrow = nsupr - nsupc; // Number of rows in the off-diagonal blocks
     
-    // Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc..., fsupc:jcol)
+    // Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc:jcol, fsupc:jcol)
     Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) );
-    Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc);
-    u = A.triangularView<Lower>().solve(u);
+//     Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc);
+    VectorBlock<ScalarVector> u(lusup, ufirst, nsupc);
+    u = A.triangularView<Lower>().solve(u); // Call the Eigen dense triangular solve interface
     
     // Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol)
     new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); 
-    Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nsupc); 
+//     Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nrow); 
+    VectorBlock<ScalarVector> l(lusup, ufirst+nsupc, nrow); 
     l = l - A * u; 
     
     return 0;