Checking Syntax...

11 files changed, 479 insertions, 205 deletions

diff --git a/Eigen/src/OrderingMethods/Eigen_Colamd.h b/Eigen/src/OrderingMethods/Eigen_Colamd.h
new file mode 100644
index 000000000..8caee7740
--- /dev/null
+++ b/Eigen/src/OrderingMethods/Eigen_Colamd.h
@@ -0,0 +1,5 @@
+
+#ifndef EIGEN_COLAMD_H
+#define EIGEN_COLAMD_H
+
+#endif
\ No newline at end of file
diff --git a/Eigen/src/OrderingMethods/Ordering.h b/Eigen/src/OrderingMethods/Ordering.h
new file mode 100644
index 000000000..c43c381a4
--- /dev/null
+++ b/Eigen/src/OrderingMethods/Ordering.h
@@ -0,0 +1,214 @@
+ 
+// This file is part of Eigen, a lightweight C++ template library
+// for linear algebra.
+//
+// Copyright (C) 2012  Désiré Nuentsa-Wakam <desire.nuentsa_wakam@inria.fr>
+//
+// Eigen is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 3 of the License, or (at your option) any later version.
+//
+// Alternatively, you can redistribute it and/or
+// modify it under the terms of the GNU General Public License as
+// published by the Free Software Foundation; either version 2 of
+// the License, or (at your option) any later version.
+//
+// Eigen is distributed in the hope that it will be useful, but WITHOUT ANY
+// WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+// FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License or the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU Lesser General Public
+// License and a copy of the GNU General Public License along with
+// Eigen. If not, see <http://www.gnu.org/licenses/>.
+
+#ifndef EIGEN_ORDERING_H
+#define EIGEN_ORDERING_H
+
+#include <Eigen_Colamd.h> 
+#include <Amd.h>
+
+namespace Eigen {
+template<class Derived> 
+class OrderingBase
+{
+  public:
+    typedef typename internal::traits<Derived>::MatrixType MatrixType;
+    typedef typename MatrixType::Scalar Scalar;
+    typedef typename MatrixType::Index Index;
+    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+  
+  public:
+    OrderingBase():m_isInitialized(false)
+    {
+      
+    }
+    OrderingBase(const MatrixType& mat):OrderingBase()
+    {
+      compute(mat);
+    }
+    Derived& compute(const MatrixType& mat)
+    {
+      return derived().compute(mat);
+    }
+    Derived& derived()
+    {
+      return *static_cast<Derived*>(this);
+    }
+    const Derived& derived() const
+    {
+      return *static_cast<const Derived*>(this);
+    }
+    /** 
+     * Get the permutation vector 
+     */
+    PermutationType& get_perm(const MatrixType& mat)
+    {
+      if (m_isInitialized = true) return m_P; 
+      else abort(); // FIXME Should find a smoother way to exit with error code
+    }
+    template<typename MatrixType> 
+    void at_plus_a(const MatrixType& mat);
+    
+    /** keeps off-diagonal entries; drops diagonal entries */
+    struct keep_diag {
+      inline bool operator() (const Index& row, const Index& col, const Scalar&) const
+      {
+        return row!=col;
+      }
+    };
+    
+  protected:
+    void init()
+    {
+      m_isInitialized = false;
+    }
+    PermutationType m_P; // The computed permutation 
+    mutable bool m_isInitialized; 
+    SparseMatrix<Scalar,ColMajor,Index> m_mat; // Stores the (symmetrized) matrix to permute
+}
+/**
+ * Get the symmetric pattern A^T+A from the input matrix A. 
+ * NOTE: The values should not be considered here
+ */
+template<typename MatrixType>
+void OrderingBase::at_plus_a(const MatrixType& mat)
+{
+    MatrixType C;
+    C = mat.transpose(); // NOTE: Could be  costly
+    for (int i = 0; i < C.rows(); i++) 
+    {
+        for (typename MatrixType::InnerIterator it(C, i); it; ++it)
+          it.valueRef() = 0.0;
+    }
+    m_mat = C + mat; 
+
+/** 
+ * Get the column approximate minimum degree ordering 
+ * The matrix should be in column-major format
+ */
+template<typename Scalar, typename Index>
+class COLAMDOrdering: public OrderingBase< ColamdOrdering<Scalar, Index> >
+{
+  public:
+    typedef OrderingBase< ColamdOrdering<Scalar, Index> > Base;
+    typedef SparseMatrix<Scalar,ColMajor,Index> MatrixType;  
+    
+  public:
+    COLAMDOrdering():Base() {}
+    
+    COLAMDOrdering(const MatrixType& matrix):Base()
+    {
+      compute(matrix);
+    }
+    COLAMDOrdering(const MatrixType& mat, PermutationType& perm_c):Base()
+    {
+      compute(matrix); 
+      perm_c = this.get_perm();
+    }
+    void compute(const MatrixType& mat)
+    {
+      // Test if the matrix is column major...
+          
+      int m = mat.rows();
+      int n = mat.cols();
+      int nnz = mat.nonZeros();
+      // Get the recommended value of Alen to be used by colamd
+      int Alen = colamd_recommended(nnz, m, n); 
+      // Set the default parameters
+      double knobs[COLAMD_KNOBS]; 
+      colamd_set_defaults(knobs);
+      
+      int info;
+      VectorXi p(n), A(nnz); 
+      for(int i=0; i < n; i++) p(i) = mat.outerIndexPtr()(i);
+      for(int i=0; i < nnz; i++) A(i) = mat.innerIndexPtr()(i);
+      // Call Colamd routine to compute the ordering 
+      info = colamd(m, n, Alen, A,p , knobs, stats)
+      eigen_assert( (info != FALSE)&& "COLAMD failed " );
+      
+      m_P.resize(n);
+      for (int i = 0; i < n; i++) m_P(p(i)) = i;
+      m_isInitialized = true;
+    }
+  protected:
+    using Base::m_isInitialized;
+    using Base m_P; 
+}
+
+/** 
+ * Get the approximate minimum degree ordering
+ * If the matrix is not structurally symmetric, an ordering of A^T+A is computed
+ * \tparam Scalar The type of the scalar of the matrix for which the ordering is applied
+ * \tparam  Index The type of indices of the matrix 
+ * \tparam _UpLo If the matrix is symmetric, indicates which part to use
+ */
+template <typename Scalar, typename Index, typename _UpLo>
+class AMDordering : public OrderingBase<AMDOrdering<Scalar, Index> >
+{
+  public:
+    enum { UpLo = _UpLo };
+    typedef OrderingBase< AMDOrdering<Scalar, Index> > Base;
+    typedef SparseMatrix<Scalar, ColMajor,Index> MatrixType;  
+  public:
+    AMDOrdering():Base(){}
+    AMDOrdering(const MatrixType& mat):Base()
+    {
+      compute(matrix);
+    }
+    AMDOrdering(const MatrixType& mat, PermutationType& perm_c):Base()
+    {
+      compute(matrix); 
+      perm_c = this.get_perm();
+    }
+    /** Compute the permutation vector from a column-major sparse matrix */
+    void compute(const MatrixType& mat)
+    {
+      // Compute the symmetric pattern
+      at_plus_a(mat); 
+    
+      // Call the AMD routine 
+      m_mat.prune(keep_diag());
+      internal::minimum_degree_ordering(m_mat, m_P);
+      if (m_P.size()>0) m_isInitialized = true;
+    }
+    /** Compute the permutation with a self adjoint matrix */
+    template <typename SrcType, unsigned int SrcUpLo> 
+    void compute(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat)
+    {
+      m_mat = mat;
+      
+      // Call the AMD routine 
+      m_mat.prune(keep_diag());
+      internal::minimum_degree_ordering(m_mat, m_P);
+      if (m_P.size()>0) m_isInitialized = true;
+    }
+  protected:
+    using Base::m_isInitialized;
+    using Base::m_P;
+    using Base::m_mat;
+}
+
+} // end namespace Eigen
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index e3838fbb7..a4b4fa98b 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -32,12 +32,22 @@ template <typename _MatrixType>
 class SparseLU;
 
 #include <Ordering.h>
-#include <SparseLU_Utils.h>
-#include <SuperNodalMatrix.h>
 #include <SparseLU_Structs.h>
 #include <SparseLU_Memory.h>
-#include <SparseLU_Coletree.h>
+#include <SparseLU_Utils.h>
+#include <SuperNodalMatrix.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
+ * 
+ * This class implements the supernodal LU factorization for general matrices. 
+ * 
+ * \tparam _MatrixType The type of the sparse matrix. It must be a column-major SparseMatrix<>
+ */
 template <typename _MatrixType>
 class SparseLU
 {
@@ -47,7 +57,7 @@ class SparseLU
     typedef typename MatrixType::Index Index; 
     typedef SparseMatrix<Scalar,ColMajor,Index> NCMatrix;
     typedef SuperNodalMatrix<Scalar, Index> SCMatrix; 
-    typedef GlobalLU_t<Scalar, Index> Eigen_GlobalLU_t;
+    typedef GlobalLU_t<Scalar, Index> LU_GlobalLU_t;
     typedef Matrix<Scalar,Dynamic,1> ScalarVector;
     typedef Matrix<Index,Dynamic,1> IndexVector;
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
@@ -58,18 +68,28 @@ class SparseLU
     }
     SparseLU(const MatrixType& matrix):SparseLU()
     {
-      
       compute(matrix);
     }
     
     ~SparseLU()
     {
-      
+      // Free all explicit dynamic pointers 
     }
     
     void analyzePattern (const MatrixType& matrix);
     void factorize (const MatrixType& matrix);
-    void compute (const MatrixType& matrix);
+    
+    /**
+     * Compute the symbolic and numeric factorization of the input sparse matrix.
+     * The input matrix should be in column-major storage. 
+     */
+    void compute (const MatrixType& matrix)
+    {
+      // Analyze 
+      analyzePattern(matrix); 
+      //Factorize
+      factorize(matrix);
+    }
     template<typename Rhs, typename Dest>
     bool SparseLU::_solve(const MatrixBase<Rhs> &b, MatrixBase<Dest> &dest) const
     
@@ -102,6 +122,13 @@ 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_dsnode_bmod (const Index jcol, const Index jsupno, const Index fsupc, 
+                      ScalarVector& dense, ScalarVector& tempv, LU_GlobalLu_t& Glu);
     
     // Variables 
     mutable ComputationInfo m_info;
@@ -113,14 +140,12 @@ class SparseLU
     SCMatrix m_Lstore; // The lower triangular matrix (supernodal)
     NCMatrix m_Ustore; // The upper triangular matrix
     PermutationType m_perm_c; // Column permutation 
-    PermutationType m_iperm_c; // Column permutation 
     PermutationType m_perm_r ; // Row permutation
-    PermutationType m_iperm_r ; // Inverse row permutation
     IndexVector m_etree; // Column elimination tree 
     
     ScalarVector m_work; // Scalar work vector 
     IndexVector m_iwork; //Index work vector 
-    static Eigen_GlobalLU_t m_Glu; // persistent data to facilitate multiple factors 
+    static LU_GlobalLU_t m_glu; // persistent data to facilitate multiple factors 
                                // should be defined as a class member
     // SuperLU/SparseLU options 
     bool m_symmetricmode;
@@ -135,7 +160,8 @@ class SparseLU
     int m_colblk; // The minimum column dimension for 2-D blocking to be used;
     int m_fillfactor; // The estimated fills factors for L and U, compared with A
     RealScalar m_diagpivotthresh; // Specifies the threshold used for a diagonal entry to be an acceptable pivot
-    int nnzL, nnzU; // Nonzeros in L and U factors 
+    int m_nnzL, m_nnzU; // Nonzeros in L and U factors 
+  
   private:
     // Copy constructor 
     SparseLU (SparseLU& ) {}
@@ -156,45 +182,56 @@ void SparseLU::initperfvalues()
 
 /** 
  * Compute the column permutation to minimize the fill-in (file amd.c )
+ * 
  *  - Apply this permutation to the input matrix - 
+ * 
  *  - Compute the column elimination tree on the permuted matrix (file Eigen_Coletree.h)
+ * 
  *  - Postorder the elimination tree and the column permutation (file Eigen_Coletree.h)
- *  - 
+ * 
  */
-template <typename MatrixType>
+template <typename MatrixType, typename OrderingType>
 void SparseLU::analyzePattern(const MatrixType& mat)
 {
-  // Compute the column permutation 
-  AMDordering amd(mat); 
-  m_perm_c = amd.get_perm_c(); 
+  
+  //TODO  It is possible as in SuperLU to compute row and columns scaling vectors to equilibrate the matrix mat.
+  
+  // Compute the fill-reducing ordering  
+  // TODO Currently, the only  available ordering method is AMD. 
+  
+  OrderingType ord(mat); 
+  m_perm_c = ord.get_perm(); 
+  //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;  //how is the permutation represented ???
+  m_mat = mat * m_perm_c;  //FIXME Check if this is valid, check as well how to permute only the index
     
   // Compute the column elimination tree of the permuted matrix 
   if (m_etree.size() == 0)  m_etree.resize(m_mat.cols());
-  internal::sp_coletree(m_mat, m_etree); 
+  LU_sp_coletree(m_mat, m_etree); 
     
   // In symmetric mode, do not do postorder here
-  if (m_symmetricmode ==  false) {
+  if (!m_symmetricmode) {
     IndexVector post, iwork; 
     // Post order etree
-    post = internal::TreePostorder(m_mat.cols(), m_etree); 
+    LU_TreePostorder(m_mat.cols(), m_etree, post); 
       
     // Renumber etree in postorder 
     iwork.resize(n+1);
     for (i = 0; i < n; ++i) iwork(post(i)) = post(m_etree(i));
-    m_etree = iwork; 
-      
-    // Postmultiply A*Pc by post, 
-    // i.e reorder the matrix according to the postorder of the etree
-    // FIXME Check if this is available : constructor from a vector
-    PermutationType post_perm(post); 
-    m_mat = m_mat * post_perm; 
-  
-    // Product of m_perm_c  and post
-    for (i = 0; i < n; ++i) iwork(i) = m_perm_c(post_perm.indices()(i));
-    m_perm_c = iwork; 
+    m_etree = iwork;
+    
+    // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
+    PermutationType post_perm(post);
+    //m_mat = m_mat * post_perm; // FIXME This should surely be in factorize()  
+  
+    // Composition of the two permutations
+    m_perm_c = m_perm_c * post_perm;
   } // end postordering 
+  
+  m_analysisIsok = true; 
 }
 
 /** 
@@ -217,36 +254,43 @@ template <typename MatrixType>
 void  SparseLU::factorize(const MatrixType& matrix)
 {
   
-  // Allocate storage common to the factor routines
-  int lwork = 0;
-  int info = LUMemInit(lwork); 
-  eigen_assert ( (info == 0) && "Unable to allocate memory for the factors"); 
+  eigen_assert(m_analysisIsok && "analyzePattern() should be called first"); 
+  eigen_assert((matrix.rows() == matrix.cols()) && "Only for squared matrices");
+  
+  // Apply the column permutation computed in analyzepattern()
+  m_mat = matrix * m_perm_c; 
+  m_mat.makeCompressed(); 
   
   int m = m_mat.rows();
   int n = m_mat.cols();
+  int nnz = m_mat.nonZeros();
   int maxpanel = m_panel_size * m;
+  // Allocate storage common to the factor routines
+  int lwork = 0;
+  int info = LUMemInit(m, n, nnz, m_work, m_iwork, lwork, m_fillratio, m_panel_size, m_maxsuper, m_rowblk, m_glu); 
+  if (info) 
+  {
+    std::cerr << "UNABLE TO ALLOCATE WORKING MEMORY\n\n" ;
+    m_factorizationIsOk = false;
+    return ; 
+  }
   
-  // Set up pointers for integer working arrays 
-  VectorBlock<IndexVector> segrep(m_iwork, 0, m); 
-//   Map<IndexVector> segrep(&m_iwork(0), m); //
-  
-  VectorBlock<IndexVector> parent(segrep, m, m); 
-//   Map<IndexVector> parent(&segrep(0) + m, m); //
-  
-  VectorBlock<IndexVector> xplore(parent, m, m); 
-//   Map<IndexVector> xplore(&parent(0) + m, m); //
-  
-  VectorBlock<IndexVector> repnfnz(xplore, m, maxpanel); 
-//   Map<IndexVector> repfnz(&xplore(0) + m, maxpanel); //
-  
-  VectorBlock<IndexVector> panel_lsub(repfnz, maxpanel, maxpanel) 
-//   Map<IndexVector> panel_lsub(&repfnz(0) + maxpanel, maxpanel);//
-  
-  VectorBlock<IndexVector> xprune(panel_lsub, maxpanel, n); 
-//   Map<IndexVector> xprune(&panel_lsub(0) + maxpanel, n); //
   
-  VectorBlock<IndexVector> marker(xprune, n, m * LU_NO_MARKER); 
-//   Map<IndexVector> marker(&xprune(0)+n, m * LU_NO_MARKER); //
+  // Set up pointers for integer working arrays 
+  int idx = 0; 
+  VectorBlock<IndexVector> segrep(m_iwork, idx, m); 
+  idx += m; 
+  VectorBlock<IndexVector> parent(m_iwork, idx, m);
+  idx += m;
+  VectorBlock<IndexVector> xplore(m_iwork, idx, m); 
+  idx += m; 
+  VectorBlock<IndexVector> repnfnz(m_iwork, idx, maxpanel);
+  idx += maxpanel;
+  VectorBlock<IndexVector> panel_lsub(m_iwork, idx, maxpanel)
+  idx += maxpanel; 
+  VectorBlock<IndexVector> xprune(m_iwork, idx, n);
+  idx += n; 
+  VectorBlock<IndexVector> marker(m_iwork, idx, m * LU_NO_MARKER); 
   
   repfnz.setConstant(-1); 
   panel_lsub.setConstant(-1);
@@ -259,43 +303,41 @@ void  SparseLU::factorize(const MatrixType& matrix)
   
   // Setup Permutation vectors
   // Compute the inverse of perm_c
-  PermutationType iperm_c;
-  iperm_c = m_perm_c.inverse();
+  PermutationType iperm_c (m_perm_c.inverse() );
   
   // Identify initial relaxed snodes
   IndexVector relax_end(n);
   if ( m_symmetricmode = true ) 
-    LU_heap_relax_snode(n, m_etree, m_relax, marker, relax_end);
+    internal::LU_heap_relax_snode(n, m_etree, m_relax, marker, relax_end);
   else
-    LU_relax_snode(n, m_etree, m_relax, marker, relax_end);
+    internal::LU_relax_snode(n, m_etree, m_relax, marker, relax_end);
   
   m_perm_r.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;
-  Index& nzlumax = m_Glu.nzlumax; 
+  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;
+  Index& nzlumax = m_glu.nzlumax; 
     
   supno(0) = IND_EMPTY; 
-  xsup(0) = xlsub(0) = xusub(0) = xlusup(0);
+  xsup(0) = xlsub(0) = xusub(0) = xlusup(0) = 0;
   int panel_size = m_panel_size; 
-  int wdef = panel_size; // upper bound on panel width
+  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
   //  (b) panel_size contiguous columns, <panel_size> defined by the user
   register int jcol,kcol; 
-  int min_mn = std::min(m,n);
   IndexVector panel_histo(n);
   Index nextu, nextlu, jsupno, fsupc, new_next;
-  int pivrow; // Pivotal row number in the original row matrix
+  Index pivrow; // Pivotal row number in the original row matrix
   int nseg1; // Number of segments in U-column above panel row jcol
   int nseg; // Number of segments in each U-column 
   int irep,ir; 
-  for (jcol = 0; jcol < min_mn; )
+  for (jcol = 0; jcol < n; )
   {
     if (relax_end(jcol) != IND_EMPTY) 
     { // Starting a relaxed node from jcol
@@ -308,6 +350,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
       {
         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
@@ -315,16 +358,17 @@ void  SparseLU::factorize(const MatrixType& matrix)
       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);
-      nzlumax = m_Glu.nzlumax;
       while (new_next > nzlumax ) 
       {
-        mem = LUMemXpand<Scalar>(lusup, nzlumax, nextlu, LUSUP, m_Glu);
+        mem = LUMemXpand<Scalar>(lusup, nzlumax, nextlu, LUSUP, m_glu);
         if (mem) 
         {
+          std::cerr << "MEMORY ALLOCATION FAILED FOR L FACTOR \n"; 
           m_factorizationIsOk = false; 
           return; 
         }
       }
+      
       // Now, left-looking factorize each column within the snode
       for (icol = jcol; icol<=kcol; icol++){
         xusub(icol+1) = nextu;
@@ -336,7 +380,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         LU_snode_bmod(icol, jsupno, fsupc, dense, tempv); 
         
         // Eliminate the current column 
-        info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r, m_iperm_c, pivrow, m_Glu); 
+        info = LU_pivotL(icol, m_diagpivotthresh, m_perm_r, m_iperm_c, pivrow, m_glu); 
         if ( info ) 
         {
           m_info = NumericalIssue; 
@@ -351,7 +395,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
       
       // Adjust panel size so that a panel won't overlap with the next relaxed snode. 
       panel_size = w_def;
-      for (k = jcol + 1; k < std::min(jcol+panel_size, min_mn); k++)
+      for (k = jcol + 1; k < std::min(jcol+panel_size, n); k++)
       {
         if (relax_end(k) != IND_EMPTY) 
         {
@@ -359,14 +403,14 @@ void  SparseLU::factorize(const MatrixType& matrix)
           break; 
         }
       }
-      if (k == min_mn) 
-        panel_size = min_mn - jcol; 
+      if (k == n) 
+        panel_size = n - jcol; 
         
       // Symbolic outer factorization on a panel of columns 
-      LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r, nseg1, dense, panel_lsub, segrep, repfnz, xprune, marker, parent, xplore, m_Glu); 
+      LU_panel_dfs(m, panel_size, jcol, m_mat, m_perm_r, 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_Glu); 
+      LU_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, j< jcol + panel_size; jj++) 
@@ -377,7 +421,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         //Depth-first-search for the current column
         VectorBlock<IndexVector> panel_lsubk(panel_lsub, k, m); //FIXME
         VectorBlock<IndexVector> repfnz_k(repfnz, k, m); //FIXME 
-        info = LU_column_dfs(m, jj, perm_r, nseg, panel_lsub(k), segrep, repfnz_k, xprune, marker, parent, xplore, m_Glu); 
+        info = LU_column_dfs(m, jj, perm_r, nseg, panel_lsub(k), segrep, repfnz_k, xprune, marker, parent, xplore, m_glu); 
         if ( !info ) 
         {
           m_info = NumericalIssue; 
@@ -387,7 +431,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         // Numeric updates to this column 
         VectorBlock<IndexVector> dense_k(dense, k, m); //FIXME 
         VectorBlock<IndexVector> segrep_k(segrep, nseg1, m) // FIXME Check the length
-        info = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_Glu); 
+        info = LU_column_bmod(jj, (nseg - nseg1), dense_k, tempv, segrep_k, repfnz_k, jcol, m_glu); 
         if ( info ) 
         {
           m_info = NumericalIssue; 
@@ -397,7 +441,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         
         // 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); 
+        info = LU_copy_to_col(jj, nseg, segrep, repfnz_k, perm_r, dense_k, m_glu); 
         if ( info ) 
         {
           m_info = NumericalIssue; 
@@ -406,7 +450,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         }
         
         // Form the L-segment 
-        info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r, iperm_c, pivrow, m_Glu);
+        info = LU_pivotL(jj, m_diagpivotthresh, m_perm_r, iperm_c, pivrow, m_glu);
         if ( info ) 
         {
           m_info = NumericalIssue; 
@@ -415,7 +459,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
         }
         
         // Prune columns (0:jj-1) using column jj
-        LU_pruneL(jj, m_perm_r, pivrow, nseg, segrep, repfnz_k, xprune, m_Glu); 
+        LU_pruneL(jj, m_perm_r, pivrow, nseg, segrep, repfnz_k, xprune, m_glu); 
         
         // Reset repfnz for this column 
         for (i = 0; i < nseg; i++)
@@ -442,17 +486,17 @@ void  SparseLU::factorize(const MatrixType& matrix)
     }
   }
   // Count the number of nonzeros in factors 
-  LU_countnz(min_mn, xprune, m_nnzL, m_nnzU, m_Glu); 
+  LU_countnz(n, xprune, m_nnzL, m_nnzU, m_glu); 
   // Apply permutation  to the L subscripts 
-  LU_fixupL(min_mn, m_perm_r, m_Glu); 
+  LU_fixupL(n, m_perm_r, m_glu); 
   
   // Free work space iwork and work 
   //...
   
   // Create supernode matrix L 
-  m_Lstore.setInfos(m, min_mn, nnzL, Glu.lusup, Glu.xlusup, Glu.lsub, Glu.xlsub, Glu.supno; Glu.xsup); 
+  m_Lstore.setInfos(m, n, m_nnzL, Glu.lusup, Glu.xlusup, Glu.lsub, Glu.xlsub, Glu.supno; Glu.xsup); 
   // Create the column major upper sparse matrix  U
-  new (&m_Ustore) Map<SparseMatrix<Scalar, ColumnMajor> > ( m, min_mn, nnzU, Glu.xusub.data(), Glu.usub.data(), Glu.ucol.data() ); //FIXME 
+  new (&m_Ustore) Map<SparseMatrix<Scalar, ColumnMajor> > ( m, n, m_nnzU, Glu.xusub.data(), Glu.usub.data(), Glu.ucol.data() ); //FIXME 
   this.m_Ustore = m_Ustore; 
   
   m_info = Success;
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
index d57048883..4c42387be 100644
--- a/Eigen/src/SparseLU/SparseLU_Coletree.h
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -49,26 +49,26 @@
   * NOTE : The matrix is supposed to be in column-major format. 
   * 
   */
-template<typename MatrixType>
-int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
+template<typename MatrixType, typename IndexVector>
+int SparseLU::LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
 {
   int nc = mat.cols(); // Number of columns 
   int nr = mat.rows(); // Number of rows 
   
-  VectorXi root(nc); // root of subtree of etree 
+  IndexVector root(nc); // root of subtree of etree 
   root.setZero();
-  VectorXi pp(nc); // disjoint sets 
+  IndexVector pp(nc); // disjoint sets 
   pp.setZero(); // Initialize disjoint sets 
-  VectorXi firstcol(nr); // First nonzero column in each row 
+  IndexVector firstcol(nr); // First nonzero column in each row 
   firstcol.setZero(); 
   
-  //Compute firstcol[row]
+  //Compute first nonzero column in each row 
   int row,col; 
   firstcol.setConstant(nc);  //for (row = 0; row < nr; firstcol(row++) = nc); 
   for (col = 0; col < nc; col++)
   {
     for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
-    { // Is it necessary to brows the whole matrix, the lower part should do the job ??
+    { // Is it necessary to browse the whole matrix, the lower part should do the job ??
       row = it.row();
       firstcol(row) = std::min(firstcol(row), col);
     }
@@ -80,7 +80,7 @@ int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
   int rset, cset, rroot; 
   for (col = 0; col < nc; col++) 
   {
-    pp(col) = cset = col; // Initially, each element is in its own set 
+    cset = pp(col) = col; // Initially, each element is in its own set //FIXME
     root(cset) = col; 
     parent(col) = nc; 
     for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
@@ -92,7 +92,7 @@ int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
       if (rroot != col) 
       {
         parent(rroot) = col; 
-        pp(cset) = cset = rset; // Get the union of cset and rset 
+        cset = pp(cset) = rset; // Get the union of cset and rset  //FIXME
         root(cset) = col; 
       }
     }
@@ -101,7 +101,8 @@ int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
 }
 
 /** Find the root of the tree/set containing the vertex i : Use Path halving */ 
-int etree_find (int i, VectorXi& pp)
+template<typename IndexVector>
+int etree_find (int i, IndexVector& pp)
 {
   int p = pp(i); // Parent 
   int gp = pp(p); // Grand parent 
@@ -116,12 +117,14 @@ int etree_find (int i, VectorXi& pp)
 }
 
 /**
-  * Post order a tree
+  * Post order a tree 
+  * \param parent Input tree
+  * \param post postordered tree
   */
-VectorXi TreePostorder(int n, VectorXi& parent)
+template<typename IndexVector>
+void SparseLU::LU_TreePostorder(int n, IndexVector& parent, IndexVector& post)
 {
-  VectorXi first_kid, next_kid; // Linked list of children 
-  VectorXi post; // postordered etree
+  IndexVector first_kid, next_kid; // Linked list of children 
   int postnum; 
   // Allocate storage for working arrays and results 
   first_kid.resize(n+1); 
@@ -140,14 +143,15 @@ VectorXi TreePostorder(int n, VectorXi& parent)
   
   // Depth-first search from dummy root vertex #n
   postnum = 0; 
-  internal::nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
+  internal::LU_nr_etdfs(n, parent, first_kid, next_kid, post, postnum);
   return post; 
 }
 /** 
   * Depth-first search from vertex n.  No recursion.
   * This routine was contributed by Cédric Doucet, CEDRAT Group, Meylan, France.
 */
-void nr_etdfs (int n, int *parent, int* first_kid, int *next_kid, int *post, int postnum)
+template<typename IndexVector>
+void 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) 
@@ -155,7 +159,7 @@ void nr_etdfs (int n, int *parent, int* first_kid, int *next_kid, int *post, int
     // No kid for the current node
     first = first_kid(current);
     
-    // no first kid for the current node
+    // no kid for the current node
     if (first == -1) 
     {
       // Numbering this node because it has no kid 
@@ -169,11 +173,12 @@ void nr_etdfs (int n, int *parent, int* first_kid, int *next_kid, int *post, int
         current = parent(current); 
         // numbering the parent node 
         post(current) = postnum++;
+        
         // Get the next kid 
         next = next_kid(current); 
       }
       // stopping criterion 
-      if (postnum==n+1) return; 
+      if (postnum == n+1) return; 
       
       // Updating current node 
       current = next; 
diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h
index a981b5436..b2888e9a0 100644
--- a/Eigen/src/SparseLU/SparseLU_Memory.h
+++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -46,43 +46,48 @@
 #ifndef EIGEN_SPARSELU_MEMORY
 #define EIGEN_SPARSELU_MEMORY
 
+#define LU_NO_MARKER 3
+#define LU_NUM_TEMPV(m,w,t,b) (std::max(m, (t+b)*w)  )
+#define IND_EMPTY (-1)
+
 #define LU_Reduce(alpha) ((alpha + 1) / 2) // i.e (alpha-1)/2 + 1
 #define LU_GluIntArray(n) (5* (n) + 5)
 #define LU_TempSpace(m, w) ( (2*w + 4 + LU_NO_MARKER) * m * sizeof(Index) \
-                                  + (w + 1) * m * sizeof(Scalar)
+                                  + (w + 1) * m * sizeof(Scalar) )
+                             
 namespace internal {
   
 /**
- * \brief  Allocate various working space failed in the numerical factorization phase.
+ * \brief  Allocate various working space for the numerical factorization phase.
  * \param m number of rows of the input matrix 
  * \param n number of columns 
  * \param annz number of initial nonzeros in the matrix 
  * \param work scalar working space needed by all factor routines
  * \param iwork Integer working space 
  * \param lwork  if lwork=-1, this routine returns an estimated size of the required memory
- * \param Glu persistent data to facilitate multiple factors : will be deleted later ??
+ * \param glu persistent data to facilitate multiple factors : will be deleted later ??
  * \return an estimated size of the required memory if lwork = -1; otherwise, return the size of actually allocated when memory allocation failed 
- * NOTE Unlike SuperLU, this routine does not allow the user to provide its own user space
+ * NOTE Unlike SuperLU, this routine does not support successive factorization with the same pattern and the row permutation
  */
 template <typename ScalarVector,typename IndexVector>
-int SparseLU::LUMemInit(int m, int n, int annz, ScalarVector& work, IndexVector& iwork, int lwork, int fillratio, GlobalLU_t& Glu)
+int LUMemInit(int m, int n, int annz, ScalarVector& work, IndexVector& iwork, int lwork, int fillratio, int panel_size, int maxsuper, int rowblk,  GlobalLU_t<ScalarVector, IndexVector>& glu)
 {
   typedef typename ScalarVector::Scalar; 
   typedef typename IndexVector::Index; 
   
-  int& num_expansions = Glu.num_expansions; //No memory expansions so far
+  int& num_expansions = glu.num_expansions; //No memory expansions so far
   num_expansions = 0; 
   // Guess the size for L\U factors 
-  Index& nzlmax = Glu.nzlmax; 
-  Index& nzumax = Glu.nzumax; 
-  Index& nzlumax = Glu.nzlumax;
+  Index& nzlmax = glu.nzlmax; 
+  Index& nzumax = glu.nzumax; 
+  Index& nzlumax = glu.nzlumax;
   nzumax = nzlumax = fillratio * annz; // estimated number of nonzeros in U 
   nzlmax  = std::max(1, m_fill_ratio/4.) * annz; // estimated  nnz in L factor
 
   // Return the estimated size to the user if necessary
-  int estimated_size;
   if (lwork == IND_EMPTY) 
   {
+    int estimated_size;
     estimated_size = LU_GluIntArray(n) * sizeof(Index)  + LU_TempSpace(m, m_panel_size)
                     + (nzlmax + nzumax) * sizeof(Index) + (nzlumax+nzumax) *  sizeof(Scalar) + n); 
     return estimated_size;
@@ -91,32 +96,33 @@ int SparseLU::LUMemInit(int m, int n, int annz, ScalarVector& work, IndexVector&
   // Setup the required space 
   
   // First allocate Integer pointers for L\U factors
-  Glu.supno.resize(n+1);
-  Glu.xlsub.resize(n+1);
-  Glu.xlusup.resize(n+1);
-  Glu.xusub.resize(n+1);
+  glu.xsup.resize(n+1);
+  glu.supno.resize(n+1);
+  glu.xlsub.resize(n+1);
+  glu.xlusup.resize(n+1);
+  glu.xusub.resize(n+1);
 
   // Reserve memory for L/U factors
-  expand<ScalarVector>(Glu.lusup, nzlumax, 0, 0, num_expansions); 
-  expand<ScalarVector>(Glu.ucol,nzumax, 0, 0, num_expansions); 
-  expand<IndexVector>(Glu.lsub,nzlmax, 0, 0, num_expansions); 
-  expand<IndexVector>(Glu.usub,nzumax, 0, 1, num_expansions);
+  expand<ScalarVector>(glu.lusup, nzlumax, 0, 0, num_expansions); 
+  expand<ScalarVector>(glu.ucol,nzumax, 0, 0, num_expansions); 
+  expand<IndexVector>(glu.lsub,nzlmax, 0, 0, num_expansions); 
+  expand<IndexVector>(glu.usub,nzumax, 0, 1, num_expansions);
   
   // Check if the memory is correctly allocated, 
-  // Should be a try... catch section here 
-  while ( !Glu.lusup.size() || !Glu.ucol.size() || !Glu.lsub.size() || !Glu.usub.size())
+  // FIXME Should be a try... catch section here 
+  while ( !glu.lusup.size() || !glu.ucol.size() || !glu.lsub.size() || !glu.usub.size())
   {
-    //otherwise reduce the estimated size and retry
+    //Reduce the estimated size and retry
     nzlumax /= 2;
     nzumax /= 2;
     nzlmax /= 2;
-    //FIXME Should be an exception here
+    
     if (nzlumax < annz ) return nzlumax; 
     
-    expand<ScalarVector>(Glu.lsup, nzlumax, 0, 0, Glu); 
-    expand<ScalarVector>(Glu.ucol, nzumax, 0, 0, Glu); 
-    expand<IndexVector>(Glu.lsub, nzlmax, 0, 0, Glu); 
-    expand<IndexVector>(Glu.usub, nzumax, 0, 1, Glu); 
+    expand<ScalarVector>(glu.lsup, nzlumax, 0, 0, num_expansions); 
+    expand<ScalarVector>(glu.ucol, nzumax, 0, 0, num_expansions); 
+    expand<IndexVector>(glu.lsub, nzlmax, 0, 0, num_expansions); 
+    expand<IndexVector>(glu.usub, nzumax, 0, 1, num_expansions); 
   }
   
   // LUWorkInit : Now, allocate known working storage
@@ -194,14 +200,14 @@ int  SparseLU::expand(VectorType& vec, int& length, int len_to_copy, bool keep_p
  * \param vec vector to expand 
  * \param [in,out]maxlen On input, previous size of vec (Number of elements to copy ). on output, new size
  * \param next current number of elements in the vector.
- * \param Glu Global data structure 
+ * \param glu Global data structure 
  * \return 0 on success, > 0 size of the memory allocated so far
  */
 template <typename IndexVector>
-int SparseLU::LUMemXpand(VectorType& vec, int& maxlen, int next, LU_MemType memtype, LU_GlobalLu_t& Glu)
+int SparseLU::LUMemXpand(VectorType& vec, int& maxlen, int next, LU_MemType memtype, LU_GlobalLu_t& glu)
 {
   int failed_size; 
-  int& num_expansions = Glu.num_expansions; 
+  int& num_expansions = glu.num_expansions; 
   if (memtype == USUB)
      failed_size = expand<IndexVector>(vec, maxlen, next, 1, num_expansions);
   else
@@ -211,19 +217,19 @@ int SparseLU::LUMemXpand(VectorType& vec, int& maxlen, int next, LU_MemType memt
     return faileld_size; 
   
   // The following code  is not really needed since maxlen is passed by reference 
-  // and correspond to the appropriate field in Glu
+  // and correspond to the appropriate field in glu
 //   switch ( mem_type ) {
 //     case LUSUP:
-//       Glu.nzlumax = maxlen;
+//       glu.nzlumax = maxlen;
 //       break; 
 //     case UCOL:
-//       Glu.nzumax = maxlen; 
+//       glu.nzumax = maxlen; 
 //       break; 
 //     case LSUB:
-//       Glu.nzlmax = maxlen; 
+//       glu.nzlmax = maxlen; 
 //       break;
 //     case USUB:
-//       Glu.nzumax = maxlen; 
+//       glu.nzumax = maxlen; 
 //       break; 
 //   }
   
diff --git a/Eigen/src/SparseLU/SparseLU_Structs.h b/Eigen/src/SparseLU/SparseLU_Structs.h
index 48fde1ada..1394eccdf 100644
--- a/Eigen/src/SparseLU/SparseLU_Structs.h
+++ b/Eigen/src/SparseLU/SparseLU_Structs.h
@@ -93,33 +93,24 @@ typedef enum {DOFACT, SamePattern, Factored} fact_t;
 typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} MemType; 
 
 
-/* Obsolete, headers for dynamically managed memory 
- \tparam VectorType can be int, real scalar or complex scalar*/
-template <typename VectorType> 
-struct ExpHeader {
-  int size; // Length of the memory that has been used */
-  VectorType *mem; // Save the current pointer of the newly allocated memory
-} ExpHeader; 
-
 template <typename ScalarVector, typename IndexVector>
 struct {
-  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 
-  IndexVector* lsub; // Compressed row indices of L rectangular supernodes. 
-  IndexVector* xlusup; // pointers to the beginning of each column in lusup
-  IndexVector* xlsub; // pointers to the beginning of each column in lsub
+  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 
+  IndexVector lsub; // Compressed row indices of L rectangular supernodes. 
+  IndexVector xlusup; // pointers to the beginning of each column in lusup
+  IndexVector xlsub; // pointers to the beginning of each column in lsub
   Index   nzlmax; // Current max size of lsub
   Index   nzlumax; // Current max size of lusup
   
-  ScalarVector*  ucol; // nonzero values of U ordered by columns 
-  IndexVector* usub; // row indices of U columns in ucol
-  IndexVector* xusub; // Pointers to the beginning of each column of U in ucol 
+  ScalarVector  ucol; // nonzero values of U ordered by columns 
+  IndexVector usub; // row indices of U columns in ucol
+  IndexVector xusub; // Pointers to the beginning of each column of U in ucol 
   Index   nzumax; // Current max size of ucol
   Index   n; // Number of columns in the matrix
   
   int   num_expansions; 
-  ExpHeader *expanders; // Deprecated... Array of pointers to 4 types of memory
 } GlobalLU_t;
 
 }// End namespace Eigen 
diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h
index 8d3d5efee..5c12b6243 100644
--- a/Eigen/src/SparseLU/SparseLU_Utils.h
+++ b/Eigen/src/SparseLU/SparseLU_Utils.h
@@ -25,15 +25,13 @@
 #ifdef EIGEN_SPARSELU_UTILS_H
 #define EIGEN_SPARSELU_UTILS_H
 
-// Number of marker arrays used in the symbolic factorization  each of size n 
-#define LU_NO_MARKER 3
-#define LU_NUM_TEMPV(m,w,t,b) (std::max(m, (t+b)*w)  )
-#define IND_EMPTY (-1)
+// Number of marker arrays used in the factorization  each of size n 
 
-void SparseLU::LU_countnz(const int n, VectorXi& xprune, int& nnzL, int& nnzU, GlobalLU_t& Glu)
+template <typename IndexVector>
+void SparseLU::LU_countnz(const int n, IndexVector& xprune, int& nnzL, int& nnzU, GlobalLU_t& Glu)
 {
- VectorXi& xsup = Glu.xsup; 
- VectorXi& xlsub = Glu.xlsub; 
+ IndexVector& xsup = Glu.xsup; 
+ IndexVector& xlsub = Glu.xlsub; 
  nnzL = 0; 
  nnzU = (Glu.xusub)(n); 
  int nnzL0 = 0; 
@@ -65,12 +63,13 @@ void SparseLU::LU_countnz(const int n, VectorXi& xprune, int& nnzL, int& nnzU, G
  * and applies permutation to the remaining subscripts
  * 
  */
-void SparseLU::LU_fixupL(const int n, const VectorXi& perm_r, GlobalLU_t& Glu)
+template <typename IndexVector>
+void SparseLU::LU_fixupL(const int n, const IndexVector& perm_r, GlobalLU_t& Glu)
 {
   int nsuper, fsupc, i, j, k, jstart; 
-  VectorXi& xsup = GLu.xsup; 
-  VectorXi& lsub = Glu.lsub; 
-  VectorXi& xlsub = Glu.xlsub; 
+  IndexVector& xsup = GLu.xsup; 
+  IndexVector& lsub = Glu.lsub; 
+  IndexVector& xlsub = Glu.xlsub; 
   
   int nextl = 0; 
   int nsuper = (Glu.supno)(n); 
diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
index 908f4d4cb..4190e0462 100644
--- a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
@@ -40,9 +40,10 @@
  * the code was modified is included with the above copyright notice.
  */
 
-#ifndef EIGEN_HEAP_RELAX_SNODE_H
-#define EIGEN_HEAP_RELAX_SNODE_H
-#include <Eigen_coletree.h>
+#ifndef SPARSELU_HEAP_RELAX_SNODE_H
+#define SPARSELU_HEAP_RELAX_SNODE_H
+#include <SparseLU_coletree.h>
+namespace internal {
 /** 
  * \brief Identify the initial relaxed supernodes
  * 
@@ -53,19 +54,20 @@
  * \param descendants Number of descendants of each node in the etree
  * \param relax_end last column in a supernode
  */
-void internal::LU_heap_relax_snode (const int n, VectorXi& et, const int relax_columns, VectorXi& descendants, VectorXi& relax_end)
+template <typename IndexVector>
+void 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  
-  // Post order etree
-  VectorXi post = internal::TreePostorder(n, et); 
-  VectorXi inv_post(n+1); 
+  IndexVector post;
+  TreePostorder(n, et, post); // Post order etree
+  IndexVector inv_post(n+1); 
   register int i;
-  for (i = 0; i < n+1; ++i) inv_post(post(i)) = i;
+  for (i = 0; i < n+1; ++i) inv_post(post(i)) = i; // inv_post = post.inverse()???
   
   // Renumber etree in postorder 
-  VectorXi iwork(n);
-  VectorXi et_save(n+1);
+  IndexVector iwork(n);
+  IndexVector et_save(n+1);
   for (i = 0; i < n; ++i)
   {
     iwork(post(i)) = post(et(i));
@@ -74,7 +76,7 @@ void internal::LU_heap_relax_snode (const int n, VectorXi& et, const int relax_c
   et = iwork; 
   
   // compute the number of descendants of each node in the etree
-  relax_end.setConstant(-1);
+  relax_end.setConstant(IND_EMPTY);
   register int j, parent; 
   descendants.setZero();
   for (j = 0; j < n; j++) 
@@ -130,4 +132,5 @@ void internal::LU_heap_relax_snode (const int n, VectorXi& et, const int relax_c
   // Recover the original etree
   et = et_save; 
 }
+} // end namespace internal
 #endif
diff --git a/Eigen/src/SparseLU/SparseLU_relax_snode.h b/Eigen/src/SparseLU/SparseLU_relax_snode.h
index 61b8e74bb..f7b478560 100644
--- a/Eigen/src/SparseLU/SparseLU_relax_snode.h
+++ b/Eigen/src/SparseLU/SparseLU_relax_snode.h
@@ -40,25 +40,26 @@
  * the code was modified is included with the above copyright notice.
  */
 
-#ifndef EIGEN_HEAP_RELAX_SNODE_H
-#define EIGEN_HEAP_RELAX_SNODE_H
-#include <Eigen_coletree.h>
+#ifndef SPARSELU_RELAX_SNODE_H
+#define SPARSELU_RELAX_SNODE_H
+namespace internal {
 /** 
  * \brief Identify the initial relaxed supernodes
  * 
- * This routine applied to a column elimination tree. 
+ * This routine is applied to a column elimination tree. 
  * It assumes that the matrix has been reordered according to the postorder of the etree
  * \param et elimination tree 
  * \param relax_columns Maximum number of columns allowed in a relaxed snode 
  * \param descendants Number of descendants of each node in the etree
  * \param relax_end last column in a supernode
  */
-void internal::LU_relax_snode (const int n, VectorXi& et, const int relax_columns, VectorXi& descendants, VectorXi& relax_end)
+template <typename IndexVector>
+void 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
   register int j, parent; 
-  relax_end.setConstant(-1);
+  relax_end.setConstant(IND_EMPTY);
   descendants.setZero();
   for (j = 0; j < n; j++) 
   {
@@ -86,4 +87,6 @@ void internal::LU_relax_snode (const int n, VectorXi& et, const int relax_column
   } // End postorder traversal of the etree
   
 }
+
+} // end namespace internal
 #endif
diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
index 9da986497..6130a5622 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_bmod.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -42,16 +42,18 @@
  * granted, provided the above notices are retained, and a notice that
  * the code was modified is included with the above copyright notice.
  */
+namespace internal {
 #ifndef SPARSELU_SNODE_BMOD_H
 #define SPARSELU_SNODE_BMOD_H
-template <typename VectorType>
-int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc, 
-                              VectorType& dense, VectorType& tempv, LU_GlobalLu_t& Glu)
+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)
 {
-  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(*)
-  VectorType& lusup = Glu.lusup; // Numerical values of the rectangular supernodes
-  VectorXi& xlusup = Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
+  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(*)
   
   int nextlu = xlusup(jcol); // Starting location of the next column to add 
   int irow, isub; 
@@ -85,4 +87,5 @@ int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc,
     
     return 0;
 }
+} // End namespace internal 
 #endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
index cf64eb747..669f172f5 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
@@ -42,8 +42,9 @@
  * granted, provided the above notices are retained, and a notice that
  * the code was modified is included with the above copyright notice.
  */
-#ifdef EIGEN_SNODE_DFS_H
-#define EIGEN_SNODE_DFS_H
+#ifdef SPARSELU_SNODE_DFS_H
+#define SPARSELU_SNODE_DFS_H
+namespace eigen {
   /**
    * \brief Determine the union of the row structures of those columns within the relaxed snode.
  *    NOTE: The relaxed snodes are leaves of the supernodal etree, therefore, 
@@ -57,15 +58,15 @@
  * \param marker (in/out) working vector
  * \return 0 on success, > 0 size of the memory when memory allocation failed
  */
-  template <typename IndexVector>
-  int SparseLU::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 IndexVector, typename ScalarVector>
+  int SparseLU::LU_snode_dfs(const int jcol, const int kcol, const IndexVector* asub, const IndexVector* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t& glu)
   {
     typedef typename IndexVector::Index; 
-    IndexVector& xsup = Glu.xsup; 
-    IndexVector& supno = Glu.supno; // Supernode number corresponding to this column
-    IndexVector& lsub = Glu.lsub;
-    IndexVector& xlsub = Glu.xlsub;
-    Index& nzlmax = Glu.nzlmax; 
+    IndexVector& xsup = glu.xsup; 
+    IndexVector& supno = glu.supno; // Supernode number corresponding to this column
+    IndexVector& lsub = glu.lsub;
+    IndexVector& xlsub = glu.xlsub;
+    Index& nzlmax = glu.nzlmax; 
     int mem; 
     Index nsuper = ++supno(jcol); // Next available supernode number
     register int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
@@ -85,7 +86,7 @@
           lsub(nextl++) = krow; 
           if( nextl >= nzlmax )
           {
-            mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, Glu);
+            mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu);
             if (mem) return mem; 
           }
         }
@@ -99,13 +100,13 @@
       Index new_next = nextl + (nextl - xlsub(jcol));
       while (new_next > nzlmax)
       {
-        mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, Glu);
+        mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu);
         if (mem) return mem; 
       }
       Index ifrom, ito = nextl; 
       for (ifrom = xlsub(jcol); ifrom < nextl;)
         lsub(ito++) = lsub(ifrom++);
-      for (i = jcol+1; i <=kcol; i++)xlsub(i) = nextl;
+      for (i = jcol+1; i <=kcol; i++) xlsub(i) = nextl;
       nextl = ito;
     }
     xsup(nsuper+1) = kcol + 1; // Start of next available supernode
@@ -115,5 +116,5 @@
     return 0;
   }
    
-  
+} // end namespace eigen
 #endif
\ No newline at end of file