Add preliminary files for SparseLU

12 files changed, 1781 insertions, 0 deletions

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
new file mode 100644
index 000000000..f5a1c787e
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -0,0 +1,341 @@
+// 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_SPARSE_LU
+#define EIGEN_SPARSE_LU
+
+#include <Ordering.h>
+#include <SparseLU_Utils.h>
+#include <SuperNodalMatrix.h>
+#include <SparseLU_Structs.h>
+#include <SparseLU_Memory.h>
+#include <SparseLU_Coletree.h>
+namespace Eigen {
+
+template <typename _MatrixType>
+class SparseLU
+{
+  public:
+    typedef _MatrixType MatrixType; 
+    typedef typename MatrixType::Scalar Scalar; 
+    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 Matrix<Scalar,Dynamic,1> VectorType;
+    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+  public:
+    SparseLU():m_isInitialized(true),m_symmetricmode(false),m_fact(DOFACT),m_diagpivotthresh(1.0)
+    {
+      initperfvalues(); 
+    }
+    SparseLU(const MatrixType& matrix):SparseLU()
+    {
+      
+      compute(matrix);
+    }
+    
+    ~SparseLU()
+    {
+      
+    }
+    
+    void analyzePattern (const MatrixType& matrix);
+    void factorize (const MatrixType& matrix);
+    void compute (const MatrixType& matrix);
+    
+    /** Indicate that the pattern of the input matrix is symmetric */
+    void isSymmetric(bool sym)
+    {
+      m_symmetricmode = sym;
+    }
+    
+    /** Set the threshold used for a diagonal entry to be an acceptable pivot. */
+    void diagPivotThresh(RealScalar thresh)
+    {
+      m_diagpivotthresh = thresh; 
+    }
+  protected:
+    // Functions 
+    void initperfvalues(); 
+    
+    // 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
+    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
+    VectorXi m_etree; // Column elimination tree 
+    
+    Scalar *m_work; //
+    Index *m_iwork; //
+    static Eigen_GlobalLU_t m_Glu; // persistent data to facilitate multiple factors 
+                               // should be defined as a class member
+    // SuperLU/SparseLU options 
+    bool m_symmetricmode;
+    
+    // values for performance 
+    int m_panel_size; // a panel consists of at most <panel_size> consecutive columns
+    int m_relax; // To control degree of relaxing supernodes. If the number of nodes (columns) 
+                 // in a subtree of the elimination tree is less than relax, this subtree is considered 
+                 // as one supernode regardless of the row structures of those columns
+    int m_maxsuper; // The maximum size for a supernode in complete LU
+    int m_rowblk; // The minimum row dimension for 2-D blocking to be used;
+    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
+    
+  private:
+    // Copy constructor 
+    SparseLU (SparseLU& ) {}
+  
+}; // End class SparseLU
+
+/* Set the  default values for performance    */
+void SparseLU::initperfvalues()
+{
+  m_panel_size = 12; 
+  m_relax = 1; 
+  m_maxsuper = 100; 
+  m_rowblk = 200; 
+  m_colblk = 60; 
+  m_fillfactor = 20; 
+}
+
+
+/** 
+ * 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>
+void SparseLU::analyzePattern(const MatrixType& mat)
+{
+  // Compute the column permutation 
+  AMDordering amd(mat); 
+  m_perm_c = amd.get_perm_c(); 
+  // Apply the permutation to the column of the input  matrix
+  m_mat = mat * m_perm_c;  //how is the permutation represented ???
+    
+  // 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); 
+    
+  // In symmetric mode, do not do postorder here
+  if (m_symmetricmode ==  false) {
+    VectorXi post, iwork; 
+    // Post order etree
+    post = internal::TreePostorder(m_mat.cols(), m_etree); 
+      
+    // 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; 
+  } // end postordering 
+}
+
+/** 
+ *  - Numerical factorization 
+ *  - Interleaved with the symbolic factorization 
+ * \tparam MatrixType The type of the matrix, it should be a column-major sparse matrix
+ * \return info where
+ *  : successful exit
+ *    = 0: successful exit
+ *    > 0: if info = i, and i is
+ *       <= A->ncol: U(i,i) is exactly zero. The factorization has
+ *          been completed, but the factor U is exactly singular,
+ *          and division by zero will occur if it is used to solve a
+ *          system of equations.
+ *       > A->ncol: number of bytes allocated when memory allocation
+ *         failure occurred, plus A->ncol. If lwork = -1, it is
+ *         the estimated amount of space needed, plus A->ncol.  
+ */
+template <typename MatrixType>
+int 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"); 
+  
+  int m = m_mat.rows();
+  int n = m_mat.cols();
+  int maxpanel = m_panel_size * m;
+  
+  // Set up pointers for integer working arrays 
+  Map<VectorXi> segrep(m_iwork, m); //
+  Map<VectorXi> parent(&segrep(0) + m, m); //
+  Map<VectorXi> xplore(&parent(0) + m, m); //
+  Map<VectorXi> repfnz(&xplore(0) + m, maxpanel); // 
+  Map<VectorXi> panel_lsub(&repfnz(0) + maxpanel, maxpanel);//
+  Map<VectorXi> xprune(&panel_lsub(0) + maxpanel, n); //
+  Map<VectorXi> marker(&xprune(0)+n, m * LU_NO_MARKER); //
+  repfnz.setConstant(-1); 
+  panel_lsub.setConstant(-1);
+  
+  // Set up pointers for scalar working arrays 
+  VectorType dense(maxpanel);
+  dense.setZero();
+  VectorType tempv(LU_NUM_TEMPV(m,m_panel_size,m_maxsuper,m_rowblk); 
+  tempv.setZero();
+  
+  // Setup Permutation vectors
+  PermutationType iperm_r; // inverse of perm_r
+  if (m_fact = SamePattern_SameRowPerm)
+    iperm_r = m_perm_r.inverse();
+  // Compute the inverse of perm_c
+  PermutationType iperm_c;
+  iperm_c = m_perm_c.inverse();
+  
+  // Identify initial relaxed snodes
+  VectorXi relax_end(n);
+  if ( m_symmetricmode = true ) 
+    LU_heap_relax_snode(n, m_etree, m_relax, marker, relax_end);
+  else
+    LU_relax_snode(n, m_etree, m_relax, marker, relax_end);
+  
+  m_perm_r.setConstant(-1);
+  marker.setConstant(-1);
+  
+  VectorXi& xsup = m_Glu.xsup; 
+  VectorXi& supno = m_GLu.supno; 
+  VectorXi& xlsub = m_Glu.xlsub;
+  VectorXi& xlusup = m_GLu.xlusup;
+  VectorXi& xusub = m_Glu.xusub;
+    
+  supno(0) = -1; 
+  xsup(0) = xlsub(0) = xusub(0) = xlusup(0);
+  int panel_size = m_panel_size; 
+  int wdef = 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);
+  VectorXi panel_histo(n);
+  bool ok = true; 
+  Index nextu, nextlu, jsupno, fsupc, new_next;
+  int 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 
+  for (jcol = 0; jcol < min_mn; )
+  {
+    if (relax_end(jcol) != -1) 
+    { // Starting a relaxed node from jcol
+      kcol = relax_end(jcol); // End index of the relaxed snode 
+      
+      // Factorize the relaxed supernode(jcol:kcol)
+      // First, determine the union of the row structure of the snode 
+      info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker); 
+      if ( !info ) 
+      {
+        ok = false; 
+        break; 
+      }
+      nextu = xusub(jcol); //starting location of column jcol in ucol
+      nextlu = xlusup(jcol); //Starting location of column jcol in lusup (rectangular supernodes)
+      jsupno = supno(jcol); // Supernode number which column jcol belongs to 
+      fsupc = xsup(jsupno); //First column number of the current supernode
+      new_next = nextlu + (xlsub(fsupc+1)-xlsub(fsupc)) * (kcol - jcol + 1);
+      nzlumax = m_Glu.nzlumax;
+      while (new_next > nzlumax ) 
+      {
+        m_Glu.lusup = LUMemXpand<Scalar>(jcol, nextlu, LUSUP, nzlumax);
+        m_GLu.nzlumax = nzlumax;
+      }
+      // Now, left-looking factorize each column within the snode
+      for (icol = jcol; icol<=kcol; icol++){
+        xusub(icol+1) = nextu;
+        // Scatter into SPA dense(*)
+        for (typename MatrixType::InnerIterator it(m_mat, icol); it; ++it)
+          dense(it.row()) = it.val();
+        
+        // Numeric update within the snode 
+        LU_snode_bmod(icol, jsupno, fsupc, dense, tempv); 
+        
+        // Eliminate the current column 
+        info = LU_pivotL(icol,  pivrow); 
+        eigen_assert(info == 0 && "The matrix is structurally singular"); 
+      }
+      jcol = icol; // The last column te be eliminated
+    }
+    else 
+    { // 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;
+      for (k = jcol + 1; k < std::min(jcol+panel_size, min_mn); k++)
+      {
+        if (relax_end(k) != -1) 
+        {
+          panel_size = k - jcol; 
+          break; 
+        }
+      }
+      if (k == min_mn) 
+        panel_size = min_mn - 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); 
+      
+      // Numeric sup-panel updates in topological order 
+      LU_panel_bmod(m, panel_size, jcol); 
+      
+      // Sparse LU within the panel, and below the panel diagonal 
+      for ( jj = jcol, j< jcol + panel_size; jj++) 
+      {
+        k = (jj - jcol) * m; // Column index for w-wide arrays 
+      } // end for 
+      jcol += panel_size;  // Move to the next panel
+    } // end else 
+  } // end for -- end elimination 
+  m_info = ok ? Success : NumericalIssue;
+  m_factorizationIsOk = ok;
+}
+
+
+} // End namespace Eigen 
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_Coletree.h b/Eigen/src/SparseLU/SparseLU_Coletree.h
new file mode 100644
index 000000000..d57048883
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_Coletree.h
@@ -0,0 +1,188 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of sp_coletree.c file in SuperLU 
+ 
+ * -- SuperLU routine (version 3.1) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * August 1, 2008
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * the code was modified is included with the above copyright notice.
+ */
+#ifndef SPARSELU_COLETREE_H
+#define SPARSELU_COLETREE_H
+
+/** Compute the column elimination tree of a sparse matrix
+  * NOTE : The matrix is supposed to be in column-major format. 
+  * 
+  */
+template<typename MatrixType>
+int LU_sp_coletree(const MatrixType& mat, VectorXi& parent)
+{
+  int nc = mat.cols(); // Number of columns 
+  int nr = mat.rows(); // Number of rows 
+  
+  VectorXi root(nc); // root of subtree of etree 
+  root.setZero();
+  VectorXi pp(nc); // disjoint sets 
+  pp.setZero(); // Initialize disjoint sets 
+  VectorXi firstcol(nr); // First nonzero column in each row 
+  firstcol.setZero(); 
+  
+  //Compute firstcol[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 ??
+      row = it.row();
+      firstcol(row) = std::min(firstcol(row), col);
+    }
+  }
+  /* Compute etree by Liu's algorithm for symmetric matrices,
+          except use (firstcol[r],c) in place of an edge (r,c) of A.
+    Thus each row clique in A'*A is replaced by a star
+    centered at its first vertex, which has the same fill. */
+  int rset, cset, rroot; 
+  for (col = 0; col < nc; col++) 
+  {
+    pp(col) = cset = col; // Initially, each element is in its own set 
+    root(cset) = col; 
+    parent(col) = nc; 
+    for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
+    { //  A sequence of interleaved find and union is performed 
+      row = firstcol(it.row());
+      if (row >= col) continue; 
+      rset = internal::etree_find(row, pp); // Find the name of the set containing row
+      rroot = root(rset);
+      if (rroot != col) 
+      {
+        parent(rroot) = col; 
+        pp(cset) = cset = rset; // Get the union of cset and rset 
+        root(cset) = col; 
+      }
+    }
+  }
+  return 0;  
+}
+
+/** Find the root of the tree/set containing the vertex i : Use Path halving */ 
+int etree_find (int i, VectorXi& pp)
+{
+  int p = pp(i); // Parent 
+  int gp = pp(p); // Grand parent 
+  while (gp != p) 
+  {
+    pp(i) = gp; // Parent pointer on find path is changed to former grand parent
+    i = gp; 
+    p = pp(i);
+    gp = pp(p);
+  }
+  return p; 
+}
+
+/**
+  * Post order a tree
+  */
+VectorXi TreePostorder(int n, VectorXi& parent)
+{
+  VectorXi first_kid, next_kid; // Linked list of children 
+  VectorXi post; // postordered etree
+  int postnum; 
+  // Allocate storage for working arrays and results 
+  first_kid.resize(n+1); 
+  next_kid.setZero(n+1);
+  post.setZero(n+1);
+  
+  // Set up structure describing children
+  int v, dad; 
+  first_kid.setConstant(-1); 
+  for (v = n-1, v >= 0; v--) 
+  {
+    dad = parent(v);
+    next_kid(v) = first_kid(dad); 
+    first_kid(dad) = v; 
+  }
+  
+  // Depth-first search from dummy root vertex #n
+  postnum = 0; 
+  internal::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)
+{
+  int current = n, first, next;
+  while (postnum != n) 
+  {
+    // No kid for the current node
+    first = first_kid(current);
+    
+    // no first kid for the current node
+    if (first == -1) 
+    {
+      // Numbering this node because it has no kid 
+      post(current) = postnum++;
+      
+      // looking for the next kid 
+      next = next_kid(current); 
+      while (next == -1) 
+      {
+        // No more kids : back to the parent node
+        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; 
+      
+      // Updating current node 
+      current = next; 
+    }
+    else 
+    {
+      current = first; 
+    }
+  }
+}
+
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_Matrix.h b/Eigen/src/SparseLU/SparseLU_Matrix.h
new file mode 100644
index 000000000..c4d56ee0a
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_Matrix.h
@@ -0,0 +1,74 @@
+// 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>
+// Copyright (C) 2012 Gael Guennebaud <gael.guennebaud@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_SPARSELU_MATRIX_H
+#define EIGEN_SPARSELU_MATRIX_H
+
+/** \ingroup SparseLU_Module
+ * \brief a class to manipulate the supernodal matrices in the SparseLU factorization
+ * 
+ * This class extends the class SparseMatrix and should contain the data to easily store 
+ * and manipulate the supernodes during the factorization and solution phase of Sparse LU. 
+ * Only the lower triangular matrix has supernodes.
+ * 
+ * NOTE : This class corresponds to the SCformat structure in SuperLU
+ * 
+ */
+
+template <typename _Scalar, typename _Index>
+class SuperNodalMatrix
+{
+  public:
+    SCMatrix()
+    {
+      
+    }
+    
+    ~SCMatrix()
+    {
+      
+    }
+    operator SparseMatrix(); 
+    
+  protected:
+    Index nnz; // Number of nonzero values 
+    Index nsupper; // Index of the last supernode
+    Scalar *nzval; //array of nonzero values packed by (supernode ??) column
+    Index *nzval_colptr; //nzval_colptr[j] Stores the location in nzval[] which starts column j 
+    Index *rowind; // Array of compressed row indices of rectangular supernodes
+    Index rowind_colptr; //rowind_colptr[j] stores the location in rowind[] which starts column j
+    Index *col_to_sup; // col_to_sup[j] is the supernode number to which column j belongs
+    Index *sup_to_col; //sup_to_col[s] points to the starting column of the s-th supernode
+    // Index *nzval_colptr corresponds to m_outerIndex in SparseMatrix
+    
+  private :
+    SuperNodalMatrix(SparseMatrix& ) {}
+};
+
+SuperNodalMatrix::operator SparseMatrix() 
+{
+  
+}
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_Memory.h b/Eigen/src/SparseLU/SparseLU_Memory.h
new file mode 100644
index 000000000..6e0fc658d
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_Memory.h
@@ -0,0 +1,242 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of [s,d,c,z]memory.c files in SuperLU 
+ 
+ * -- SuperLU routine (version 3.1) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * August 1, 2008
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * the code was modified is included with the above copyright notice.
+ */
+
+#ifndef EIGEN_SPARSELU_MEMORY
+#define EIGEN_SPARSELU_MEMORY
+
+#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)
+namespace internal {
+  
+/* Allocate various working space needed in the numerical factorization phase.
+ * m_work : space fot the output data structures (lwork is the size)
+ * m_Glu: persistent data to facilitate multiple factors : is it really necessary ??
+ * NOTE Unlike SuperLU, this routine does not allow the user to provide the size to allocate 
+ * nor it return an estimated amount of space required.
+ * 
+ * Useful  variables :
+ * -  m_fillratio : Ratio of fill expected 
+ * - lwork = -1 : return an estimated size of the required memory
+ *           = 0 : Estimate and allocate the memory
+ */
+template <typename Scalar,typename Index>
+int SparseLU::LUMemInit(int lwork)
+{
+  int iword = sizeof(Index); 
+  int dword = sizeof(Scalar);
+  int n = m_Glu.n = m_mat.cols();
+  int m = m_mat.rows();
+  m_Glu.num_expansions = 0; // No memory expansions so far ??
+  int estimated_size;
+  
+  
+  if (!m_Glu.expanders)
+    m_Glu.expanders = new ExpHeader(NO_MEMTYPE); 
+  
+  if (m_fact_t != SamePattern_SameRowPerm) // Create space for a new factorization
+  { 
+    // Guess the size for L\U factors 
+    int annz = m_mat.nonZeros();
+    int nzlmax, nzumax, nzlumax; 
+    nzumax = nzlumax = m_fillratio * annz; // ???
+    nzlmax  = std::max(1, m_fill_ratio/4.) * annz; //???
+  
+    // Return the estimated size to the user if necessary
+    if (lwork = -1) 
+    {
+      estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
+                      + (nzlmax + nzumax) * iword + (nzlumax+nzumax) * dword + n); 
+      return estimated_size;
+    }
+    
+    // Setup the required space 
+    // NOTE: In SuperLU, there is an option to let the user provide its own space.
+    
+    // Allocate Integer pointers for L\U factors.resize(n+1);
+    m_Glu.supno.resize(n+1);
+    m_Glu.xlsub.resize(n+1);
+    m_Glu.xlusup.resize(n+1);
+    m_Glu.xusub.resize(n+1);
+  
+    // Reserve memory for L/U factors
+    m_Glu.lusup = internal::expand<Scalar>(nzlumax, LUSUP, 0, 0, m_Glu); 
+    m_Glu.ucol = internal::expand<Scalar>(nzumax, UCOL, 0, 0, m_Glu); 
+    m_Glu.lsub = internal::expand<Index>(nzlmax, LSUB, 0, 0, m_Glu); 
+    m_Glu.usub = internal::expand<Index>(nzumax, USUB, 0, 1, m_Glu); 
+    
+    // Check if the memory is correctly allocated, 
+    while ( !m_Glu.lusup || !m_Glu.ucol || !m_Glu.lsub || !m_Glu.usub)
+    {
+      //otherwise reduce the estimated size and retry
+      delete [] m_Glu.lusup; 
+      delete [] m_Glu.ucol;
+      delete [] m_Glu.lsub;
+      delete [] m_Glu.usub;
+      
+      nzlumax /= 2;
+      nzumax /= 2;
+      nzlmax /= 2;
+      eigen_assert (nzlumax > annz && "Not enough memory to perform factorization");
+      
+      m_Glu.lusup = internal::expand<Scalar>(nzlumax, LUSUP, 0, 0, m_Glu); 
+      m_Glu.ucol = internal::expand<Scalar>(nzumax, UCOL, 0, 0, m_Glu); 
+      m_Glu.lsub = internal::expand<Index>(nzlmax, LSUB, 0, 0, m_Glu); 
+      m_Glu.usub = internal::expand<Index>(nzumax, USUB, 0, 1, m_Glu); 
+    }
+  }
+  else  // m_fact == SamePattern_SameRowPerm;
+  {
+    if (lwork = -1) 
+    {
+      estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
+                      + (Glu.nzlmax + Glu.nzumax) * iword + (Glu.nzlumax+Glu.nzumax) * dword + n); 
+      return estimated_size;
+    }
+    //  Use existing space from previous factorization
+    // Unlike in SuperLU, this should not  be necessary here since m_Glu is persistent as a member of the class
+    m_Glu.xsup = m_Lstore.sup_to_col; 
+    m_Glu.supno = m_Lstore.col_to_sup;
+    m_Glu.xlsub = m_Lstore.rowind_colptr;
+    m_Glu.xlusup = m_Lstore.nzval_colptr;
+    xusub = m_Ustore.outerIndexPtr();
+    
+    m_Glu.expanders[LSUB].size = m_Glu.nzlmax; // Maximum value from previous factorization
+    m_Glu.expanders[LUSUP].size = m_Glu.nzlumax;
+    m_Glu.expanders[USUB].size = GLu.nzumax;
+    m_Glu.expanders[UCOL].size = m_Glu.nzumax;
+    m_Glu.lsub = GLu.expanders[LSUB].mem = m_Lstore.rowind;
+    m_Glu.lusup = GLu.expanders[LUSUP].mem = m_Lstore.nzval;
+    GLu.usub = m_Glu.expanders[USUB].mem = m_Ustore.InnerIndexPtr();
+    m_Glu.ucol = m_Glu.expanders[UCOL].mem = m_Ustore.valuePtr();
+  }
+  
+  // LUWorkInit : Now, allocate known working storage
+  int isize = (2 * m_panel_size + 3 + LU_NO_MARKER) * m + n;
+  int dsize = m * m_panel_size + LU_NUM_TEMPV(m, m_panel_size, m_maxsuper, m_rowblk); 
+  m_iwork = new Index(isize);
+  eigen_assert( (m_iwork != 0) && "Malloc fails for iwork");
+  m_work = new Scalar(dsize);
+  eigen_assert( (m_work != 0) && "Malloc fails for dwork");
+  
+  ++m_Glu.num_expansions;
+  return 0;
+} // end LuMemInit
+
+/** 
+  * Expand the existing storage to accomodate more fill-ins
+  */
+template <typename DestType >
+DestType* SparseLU::expand(int& prev_len, // Length from previous call
+        MemType type, // Which part of the memory to expand
+        int len_to_copy, // Size of the memory to be copied to new store
+        int keep_prev) // = 1: use prev_len; Do not expand this vector
+                      // = 0: compute new_len to expand)
+{
+  
+  float alpha = 1.5; // Ratio of the memory increase 
+  int new_len; // New size of the allocated memory
+  if(m_Glu.num_expansions == 0 || keep_prev) 
+    new_len = prev_len;
+  else 
+    new_len = alpha * prev_len;
+  
+  // Allocate new space 
+  DestType *new_mem, *old_mem; 
+  new_mem = new DestType(new_len); 
+  if ( m_Glu.num_expansions != 0 ) // The memory has been expanded before
+  {
+    int tries = 0; 
+    if (keep_prev) 
+    {
+      if (!new_mem) return 0; 
+    }
+    else 
+    {
+      while ( !new_mem)
+      {
+        // Reduce the size and allocate again
+        if ( ++tries > 10) return 0; 
+        alpha = LU_Reduce(alpha); 
+        new_len = alpha * prev_len; 
+        new_mem = new DestType(new_len);
+      }
+    } // keep_prev
+    //Copy the previous values to the newly allocated space 
+    ExpHeader<DestType>* expanders = m_Glu.expanders;
+    std::memcpy(new_mem, expanders[type].mem, len_to_copy); 
+    delete [] expanders[type].mem; 
+  }
+  expanders[type].mem = new_mem; 
+  expanders[type].size = new_len;
+  prev_len = new_len;
+  if(m_Glu.num_expansions) ++m_Glu.num_expansions;
+  return expanders[type].mem; 
+}
+
+/** 
+ * \brief Expand the existing storage 
+ * 
+ * NOTE: The calling sequence of this function is different from that of SuperLU
+ * 
+ * \return a pointer to the newly allocated space
+ */
+template <typename DestType>
+DestType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen)
+{
+  DestType *newmem; 
+  if (memtype == USUB)
+    new_mem = expand<DestType>(maxlen, mem_type, next, 1);
+  else
+    new_mem = expand<DestType>(maxlen, mem_type, next, 0);
+  eigen_assert(new_mem && "Can't expand memory");
+  
+  return new_mem;
+  
+}
+    
+}// Namespace Internal
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_Structs.h b/Eigen/src/SparseLU/SparseLU_Structs.h
new file mode 100644
index 000000000..72e1db343
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_Structs.h
@@ -0,0 +1,122 @@
+// 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/>.
+
+/* 
+ * NOTE: Part of this file is the modified version of files slu_[s,d,c,z]defs.h
+ * -- SuperLU routine (version 4.1) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * November, 2010
+ * 
+ * Global data structures used in LU factorization -
+ * 
+ *   nsuper: #supernodes = nsuper + 1, numbered [0, nsuper].
+ *   (xsup,supno): supno[i] is the supernode no to which i belongs;
+ *  xsup(s) points to the beginning of the s-th supernode.
+ *  e.g.   supno 0 1 2 2 3 3 3 4 4 4 4 4   (n=12)
+ *          xsup 0 1 2 4 7 12
+ *  Note: dfs will be performed on supernode rep. relative to the new 
+ *        row pivoting ordering
+ *
+ *   (xlsub,lsub): lsub[*] contains the compressed subscript of
+ *  rectangular supernodes; xlsub[j] points to the starting
+ *  location of the j-th column in lsub[*]. Note that xlsub 
+ *  is indexed by column.
+ *  Storage: original row subscripts
+ *
+ *      During the course of sparse LU factorization, we also use
+ *  (xlsub,lsub) for the purpose of symmetric pruning. For each
+ *  supernode {s,s+1,...,t=s+r} with first column s and last
+ *  column t, the subscript set
+ *    lsub[j], j=xlsub[s], .., xlsub[s+1]-1
+ *  is the structure of column s (i.e. structure of this supernode).
+ *  It is used for the storage of numerical values.
+ *  Furthermore,
+ *    lsub[j], j=xlsub[t], .., xlsub[t+1]-1
+ *  is the structure of the last column t of this supernode.
+ *  It is for the purpose of symmetric pruning. Therefore, the
+ *  structural subscripts can be rearranged without making physical
+ *  interchanges among the numerical values.
+ *
+ *  However, if the supernode has only one column, then we
+ *  only keep one set of subscripts. For any subscript interchange
+ *  performed, similar interchange must be done on the numerical
+ *  values.
+ *
+ *  The last column structures (for pruning) will be removed
+ *  after the numercial LU factorization phase.
+ *
+ *   (xlusup,lusup): lusup[*] contains the numerical values of the
+ *  rectangular supernodes; xlusup[j] points to the starting
+ *  location of the j-th column in storage vector lusup[*]
+ *  Note: xlusup is indexed by column.
+ *  Each rectangular supernode is stored by column-major
+ *  scheme, consistent with Fortran 2-dim array storage.
+ *
+ *   (xusub,ucol,usub): ucol[*] stores the numerical values of
+ *  U-columns outside the rectangular supernodes. The row
+ *  subscript of nonzero ucol[k] is stored in usub[k].
+ *  xusub[i] points to the starting location of column i in ucol.
+ *  Storage: new row subscripts; that is subscripts of PA.
+ */
+#ifndef EIGEN_LU_STRUCTS
+#define EIGEN_LU_STRUCTS
+namespace Eigen {
+
+#define NO_MEMTYPE 4 /* 0: lusup
+                        1: ucol 
+                        2: lsub 
+                        3: usub */
+typedef enum {NATURAL, MMD_ATA, MMD_AT_PLUS_A, COLAMD, MY_PREMC} colperm_t; 
+typedef enum {DOFACT, SamePattern, SamePattern_SameRowPerm, Factored} fact_t; 
+typedef enum {LUSUP, UCOL, LSUB, USUB, LLVL, ULVL} MemType; 
+
+/** Headers for dynamically managed memory 
+ \tparam BaseType can be int, real scalar or complex scalar*/
+template <typename BaseType> 
+struct ExpHeader {
+  int size; // Length of the memory that has been used */
+  BaseType *mem; 
+} ExpHeader; 
+
+template <typename VectorType, typename Index>
+struct {
+  VectorXi xsup; // supernode and column mapping 
+  VectorXi supno; // Supernode number corresponding to this column
+  VectorXi lsub; // Compressed L subscripts of rectangular supernodes
+  VectorXi xlsub; // xlsub(j) points to the starting location of the j-th column in lsub
+  VectorXi xlusup;
+  VectorXi xusub; 
+  VectorType  lusup; // L supernodes
+  VectorType  ucol; // U columns 
+  Index   nzlmax; // Current max size of lsub
+  Index   nzumax; // Current max size of ucol
+  Index   nzlumax; // Current max size of lusup
+  Index   n; // Number of columns in the matrix
+  int   num_expansions; 
+  ExpHeader *expanders; // Array of pointers to 4 types of memory
+} GlobalLU_t;
+
+}// End namespace Eigen 
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_Utils.h b/Eigen/src/SparseLU/SparseLU_Utils.h
new file mode 100644
index 000000000..3c3b24a15
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_Utils.h
@@ -0,0 +1,32 @@
+// 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/>.
+
+#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 LU_EMPTY (-1)
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
new file mode 100644
index 000000000..908f4d4cb
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_heap_relax_snode.h
@@ -0,0 +1,133 @@
+// 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/>.
+
+/* This file is a modified version of heap_relax_snode.c file in SuperLU
+ * -- SuperLU routine (version 3.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * October 15, 2003
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * 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>
+/** 
+ * \brief Identify the initial relaxed supernodes
+ * 
+ * This routine applied to a symmetric 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_heap_relax_snode (const int n, VectorXi& et, const int relax_columns, VectorXi& descendants, VectorXi& 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); 
+  register int i;
+  for (i = 0; i < n+1; ++i) inv_post(post(i)) = i;
+  
+  // Renumber etree in postorder 
+  VectorXi iwork(n);
+  VectorXi et_save(n+1);
+  for (i = 0; i < n; ++i)
+  {
+    iwork(post(i)) = post(et(i));
+  }
+  et_save = et; // Save the original etree
+  et = iwork; 
+  
+  // compute the number of descendants of each node in the etree
+  relax_end.setConstant(-1);
+  register int j, parent; 
+  descendants.setZero();
+  for (j = 0; j < n; j++) 
+  {
+    parent = et(j);
+    if (parent != n) // not the dummy root
+      descendants(parent) += descendants(j) + 1;
+  }
+  
+  // Identify the relaxed supernodes by postorder traversal of the etree
+  register int snode_start; // beginning of a snode 
+  register int k;
+  int nsuper_et_post = 0; // Number of relaxed snodes in postordered etree 
+  int nsuper_et = 0; // Number of relaxed snodes in the original etree 
+  for (j = 0; j < n; )
+  {
+    parent = et(j);
+    snode_start = j; 
+    while ( parent != n && descendants(parent) < relax_columns ) 
+    {
+      j = parent; 
+      parent = et(j);
+    }
+    // Found a supernode in postordered etree, j is the last column 
+    ++nsuper_et_post;
+    k = n;
+    for (i = snode_start; i <= j; ++i)
+      k = std::min(k, inv_post(i));
+    l = inv_post(j);
+    if ( (l - k) == (j - snode_start) )  // Same number of columns in the snode
+    {
+      // This is also a supernode in the original etree
+      relax_end(k) = l; // Record last column 
+      ++nsuper_et; 
+    }
+    else 
+    {
+      for (i = snode_start; i <= j; ++i) 
+      {
+        l = inv_post(i);
+        if (descendants(i) == 0) 
+        {
+          relax_end(l) = l;
+          ++nsuper_et;
+        }
+      }
+    }
+    j++;
+    // Search for a new leaf
+    while (descendants(j) != 0 && j < n) j++;
+  } // End postorder traversal of the etree
+  
+  // Recover the original etree
+  et = et_save; 
+}
+#endif
diff --git a/Eigen/src/SparseLU/SparseLU_panel_dfs.h b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
new file mode 100644
index 000000000..550544d05
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_panel_dfs.h
@@ -0,0 +1,221 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of xpanel_dfs.c file in SuperLU 
+ 
+ * -- SuperLU routine (version 2.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * November 15, 1997
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * the code was modified is included with the above copyright notice.
+ */
+#ifndef SPARSELU_PANEL_DFS_H
+#define SPARSELU_PANEL_DFS_H
+/**
+ * \brief Performs a symbolic factorization on a panel of columns [jcol, jcol+w)
+ * 
+ * A supernode representative is the last column of a supernode.
+ * The nonzeros in U[*,j] are segments that end at supernodes representatives
+ * 
+ * The routine returns a list of the supernodal representatives 
+ * in topological order of the dfs that generates them. This list is 
+ * a superset of the topological order of each individual column within 
+ * the panel.
+ * The location of the first nonzero in each supernodal segment 
+ * (supernodal entry location) is also returned. Each column has 
+ * a separate list for this purpose. 
+ * 
+ * Two markers arrays are used for dfs :
+ *    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 
+ * 
+ */
+template <typename MatrixType, typename VectorType>
+int SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, VectorXi& nseg, int& nseg, VectorType& dense,  VectorXi& panel_lsub, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& 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 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); 
+  nseg = 0; 
+  VectorXi& xsup = Glu.xsup; 
+  VectorXi& supno = Glu.supno; 
+  VectorXi& lsub = Glu.lsub; 
+  VectorXi& 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
+    
+    
+    // For each nnz in A[*, jj] do depth first search
+    for (MatrixType::InnerIterator it(A, jj); it; ++it)
+    {
+      krow = it.row(); 
+      dense_col(krow) = it.val(); 
+      kmark = marker(krow); 
+      if (kmark == jj) 
+        continue; // krow visited before, go to the next nonzero
+      
+      // For each unmarked krow of jj
+      marker(krow) = jj; 
+      kperm = perm_r(krow); 
+      if (kperm == -1 ) {
+        // krow is in L : place it in structure of L(*, jj)
+        panel_lsub(nextl_col++) = krow;  // krow is indexed into A
+      }
+      else 
+      {
+        // krow is in U : if its supernode-representative krep
+        // has been explored, update repfnz(*)
+        krep = xsup(supno(kperm)+1) - 1; 
+        myfnz = repfnz_col(krep); 
+        
+        if (myfnz != -1 )
+        {
+          // Representative visited before
+          if (myfnz > kperm ) repfnz_col(krep) = kperm; 
+          
+        }
+        else 
+        {
+          // Otherwise, perform dfs starting at krep
+          oldrep = -1; 
+          parent(krep) = oldrep; 
+          repfnz_col(krep) = kperm; 
+          xdfs =  xlsub(krep); 
+          maxdfs = xprune(krep); 
+          
+          do 
+          {
+            // For each unmarked kchild of krep
+            while (xdfs < maxdfs) 
+            {
+              kchild = lsub(xdfs); 
+              xdfs++; 
+              chmark = marker(kchild); 
+              
+              if (chmark != jj ) 
+              {
+                marker(kchild) = jj; 
+                chperm = perm_r(kchild); 
+                
+                if (chperm == -1) 
+                {
+                  // case kchild is in L: place it in L(*, j)
+                  panel_lsub(nextl_col++) = kchild; 
+                }
+                else
+                {
+                  // case kchild is in U :
+                  // chrep = its supernode-rep. If its rep has been explored, 
+                  // update its repfnz(*)
+                  chrep = xsup(supno(chperm)+1) - 1; 
+                  myfnz = repfnz_col(chrep); 
+                  
+                  if (myfnz != -1) 
+                  { // Visited before 
+                    if (myfnz > chperm) 
+                      repfnz_col(chrep) = chperm; 
+                  }
+                  else 
+                  { // Cont. dfs at snode-rep of kchild
+                    xplore(krep) = xdfs; 
+                    oldrep = krep; 
+                    krep = chrep; // Go deeper down G(L)
+                    parent(krep) = oldrep; 
+                    repfnz_col(krep) = chperm; 
+                    xdfs = xlsub(krep); 
+                    maxdfs = xprune(krep); 
+                    
+                  } // end if myfnz != -1
+                } // end if chperm == -1 
+                    
+              } // end if chmark !=jj
+            } // end while xdfs < maxdfs
+            
+            // krow has no more unexplored nbrs :
+            //    Place snode-rep krep in postorder DFS, if this 
+            //    segment is seen for the first time. (Note that 
+            //    "repfnz(krep)" may change later.)
+            //    Baktrack dfs to its parent
+            if (marker1(krep) < jcol )
+            {
+              segrep(nseg) = krep; 
+              ++nseg; 
+              marker1(krep) = jj; 
+            }
+            
+            kpar = parent(krep); // Pop recursion, mimic recursion 
+            if (kpar == -1) 
+              break; // dfs done 
+            krep = kpar; 
+            xdfs = xplore(krep); 
+            maxdfs = xprune(krep); 
+
+          } while (kpar != -1); // Do until empty stack 
+          
+        } // end if (myfnz = -1)
+
+      } // end if (kperm == -1) 
+
+    }// end for nonzeros in column jj
+    
+  } // end for column jj
+  
+}
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h
new file mode 100644
index 000000000..f939ef939
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
@@ -0,0 +1,132 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of dpivotL.c file in SuperLU 
+ 
+ * -- SuperLU routine (version 3.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * October 15, 2003
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * the code was modified is included with the above copyright notice.
+ */
+#ifndef SPARSELU_PIVOTL_H
+#define SPARSELU_PIVOTL_H
+/**
+ * \brief Performs the numerical pivotin on the current column of L, and the CDIV operation.
+ * 
+ * Here is the pivot policy :
+ * (1) 
+ * 
+ * \param jcol The current column of L
+ * \param pivrow [out] The pivot row
+ * 
+ * 
+ */
+int SparseLU::LU_pivotL(const int jcol, Index& pivrow)
+{
+  // Initialize pointers 
+  VectorXi& lsub = m_Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
+  VectorXi& xlsub = m_Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
+  Scalar* lusup = m_Glu.lusup.data(); // Numerical values of the rectangular supernodes
+  VectorXi& xlusup = m_Glu.xlusup; // xlusup[j] is the starting location of the j-th column in lusup(*)
+  
+  Index fsupc = (m_Glu.xsup)((m_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
+  Index* lsub_ptr = &(lsub.data()[lptr]); // Start of row indices of the supernode
+  
+  // Determine the largest abs numerical value for partial pivoting 
+  Index diagind = m_iperm_c(jcol); // diagonal index 
+  Scalar pivmax = 0.0; 
+  Index pivptr = nsupc; 
+  Index diag = -1; 
+  Index old_pivptr = nsupc; 
+  Scalar rtemp;
+  for (isub = nsupc; isub < nsupr; ++isub) {
+    rtemp = std::abs(lu_col_ptr[isub]);
+    if (rtemp > pivmax) {
+      pivmax = rtemp; 
+      pivptr = isub;
+    } 
+    if (lsub_ptr[isub] == diagind) diag = isub;
+  }
+  
+  // Test for singularity
+  if ( pivmax == 0.0 ) {
+    pivrow = lsub_ptr[pivptr];
+    m_perm_r(pivrow) = jcol;
+    return (jcol+1);
+  }
+  
+  Scalar thresh = m_diagpivotthresh * pivmax; 
+  
+  // Choose appropriate pivotal element 
+  
+  {
+    // Test if the diagonal element can be used as a pivot (given the threshold value)
+    if (diag >= 0 ) 
+    {
+      // Diagonal element exists
+      rtemp = std::abs(lu_col_ptr[diag]);
+      if (rtemp != Scalar(0.0) && rtemp >= thresh) pivptr = diag;
+    }
+    pivrow = lsub_ptr[pivptr];
+  }
+  
+  // Record pivot row
+  perm_r(pivrow) = jcol; 
+  // Interchange row subscripts
+  if (pivptr != 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++)
+    {
+      itemp = pivptr + icol * nsupr; 
+      std::swap(lu_sup_ptr[itemp], lu_sup_ptr[nsupc + icol * nsupr]);
+    }
+  }
+  // cdiv operations
+  Scalar temp = Scalar(1.0) / lu_col_ptr[nsupc];
+  for (k = nsupc+1; k < nsupr; k++)
+    lu_col_ptr[k] *= temp; 
+  return 0;
+}
+#endif
\ No newline at end of file
diff --git a/Eigen/src/SparseLU/SparseLU_relax_snode.h b/Eigen/src/SparseLU/SparseLU_relax_snode.h
new file mode 100644
index 000000000..61b8e74bb
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_relax_snode.h
@@ -0,0 +1,89 @@
+// 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/>.
+
+/* This file is a modified version of heap_relax_snode.c file in SuperLU
+ * -- SuperLU routine (version 3.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * October 15, 2003
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * 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>
+/** 
+ * \brief Identify the initial relaxed supernodes
+ * 
+ * This routine 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)
+{
+  
+  // compute the number of descendants of each node in the etree
+  register int j, parent; 
+  relax_end.setConstant(-1);
+  descendants.setZero();
+  for (j = 0; j < n; j++) 
+  {
+    parent = et(j);
+    if (parent != n) // not the dummy root
+      descendants(parent) += descendants(j) + 1;
+  }
+  
+  // Identify the relaxed supernodes by postorder traversal of the etree
+  register int snode_start; // beginning of a snode 
+  for (j = 0; j < n; )
+  {
+    parent = et(j);
+    snode_start = j; 
+    while ( parent != n && descendants(parent) < relax_columns ) 
+    {
+      j = parent; 
+      parent = et(j);
+    }
+    // Found a supernode in postordered etree, j is the last column 
+    relax_end(snode_start) = j; // Record last column
+    j++;
+    // Search for a new leaf
+    while (descendants(j) != 0 && j < n) j++;
+  } // End postorder traversal of the etree
+  
+}
+#endif
diff --git a/Eigen/src/SparseLU/SparseLU_snode_bmod.h b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
new file mode 100644
index 000000000..fc6ffc320
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_snode_bmod.h
@@ -0,0 +1,88 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of dsnode_bmod.c file in SuperLU 
+ 
+ * -- SuperLU routine (version 3.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * October 15, 2003
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * granted, provided the above notices are retained, and a notice that
+ * the code was modified is included with the above copyright notice.
+ */
+#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)
+{
+  VectorXi& lsub = m_Glu.lsub; // Compressed row subscripts of ( rectangular supernodes ??)
+  VectorXi& xlsub = m_Glu.xlsub; // xlsub[j] is the starting location of the j-th column in lsub(*)
+  Scalar* lusup = m_Glu.lusup.data(); // Numerical values of the rectangular supernodes
+  VectorXi& xlusup = m_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; 
+  // Process the supernodal portion of L\U[*,jcol]
+  for (int isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
+  {
+    irow = lsub(isub); 
+    lusup(nextlu) = dense(irow);
+    dense(irow) = 0;
+    ++nextlu;
+  }
+  xlusup(jcol + 1) = nextlu; // Initialize xlusup for next column ( jcol+1 )
+  
+  if (fsupc < jcol ){
+    int luptr = xlusup(fsupc); // points to the first column of the supernode
+    int nsupr = xlsub(fsupc + 1) -xlsub(fsupc); //Number of rows in the supernode
+    int nsupc = jcol - fsupc; // Number of columns in the supernodal portion of L\U[*,jcol]
+    int ufirst = xlusup(jcol); // points to the beginning of column jcol in supernode L\U(jsupno)
+    
+    int nrow = nsupr - nsupc; // Number of rows in the off-diagonal blocks
+    int incx = 1, incy = 1; 
+    Scalar alpha = Scalar(-1.0); 
+    Scalar beta = Scalar(1.0);
+    // Solve the triangular system for U(fsupc:jcol, jcol) with L(fspuc..., fsupc:jcol)
+    //BLASFUNC(trsv)("L", "N", "U", &nsupc, &(lusup[luptr]), &nsupr, &(lusup[ufirst]), &incx);
+    Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup[luptr]), nsupc, nsupc, OuterStride<>(nsupr) );
+    Map<Matrix<Scalar,Dynamic,1> > l(&(lusup[ufirst]), nsupc);
+    l = A.triangularView<Lower>().solve(l);
+    // Update the trailing part of the column jcol U(jcol:jcol+nrow, jcol) using L(jcol:jcol+nrow, fsupc:jcol) and U(fsupc:jcol)
+    BLASFUNC(gemv)("N", &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr, &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy);
+    
+    return 0;
+}
+#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
new file mode 100644
index 000000000..c3048be54
--- /dev/null
+++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
@@ -0,0 +1,119 @@
+// 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/>.
+
+/* 
+ 
+ * NOTE: This file is the modified version of dsnode_dfs.c file in SuperLU 
+ 
+ * -- SuperLU routine (version 2.0) --
+ * Univ. of California Berkeley, Xerox Palo Alto Research Center,
+ * and Lawrence Berkeley National Lab.
+ * November 15, 1997
+ *
+ * Copyright (c) 1994 by Xerox Corporation.  All rights reserved.
+ *
+ * THIS MATERIAL IS PROVIDED AS IS, WITH ABSOLUTELY NO WARRANTY
+ * EXPRESSED OR IMPLIED.  ANY USE IS AT YOUR OWN RISK.
+ *
+ * Permission is hereby granted to use or copy this program for any
+ * purpose, provided the above notices are retained on all copies.
+ * Permission to modify the code and to distribute modified code is
+ * 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
+  /**
+   * \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, 
+ *    the portion outside the rectangular supernode must be zero.
+ * 
+ * \param jcol start of the supernode
+ * \param kcol end of the supernode
+ * \param asub Row indices
+ * \param colptr Pointer to the beginning of each column
+ * \param xprune (out) The pruned tree ??
+ * \param marker (in/out) working vector
+ */
+  template <typename Index>
+  int SparseLU::LU_snode_dfs(const int jcol, const int kcol, const VectorXi* asub, const VectorXi* colptr, 
+                    VectorXi& xprune, VectorXi& marker, LU_GlobalLu_t *m_Glu)
+  {
+    VectorXi& xsup = m_Glu->xsup; 
+    VectorXi& supno = m_Glu->supno; // Supernode number corresponding to this column
+    VectorXi& lsub = m_Glu->lsub;
+    VectorXi& xlsub = m_Glu->xlsub;
+    
+    int nsuper = ++supno(jcol); // Next available supernode number
+    register int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
+    register int i,k; 
+    int krow,kmark; 
+    for (i = jcol; i <=kcol; i++)
+    {
+      // For each nonzero in A(*,i)
+      for (k = colptr(i); k < colptr(i+1); k++)
+      {
+        krow = asub(k); 
+        kmark = marker(krow);
+        if ( kmark != kcol )
+        {
+          // First time to visit krow
+          marker(krow) = kcol; 
+          lsub(nextl++) = krow; 
+          if( nextl >= nzlmax )
+          {
+            m_Glu->lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax);
+            m_Glu->nzlmax = nzlmax;
+            lsub = m_Glu->lsub;
+          }
+        }
+      }
+      supno(i) = nsuper; 
+    }
+    
+    // If supernode > 1, then make a copy of the subscripts for pruning
+    if (jcol < kcol)
+    {
+      int new_next = nextl + (nextl - xlsub(jcol));
+      while (new_next > nzlmax)
+      {
+        m_Glu->lsub = LUMemXpand<Index>(jcol, nextl, LSUB, &nzlmax);
+        m_Glu->nzlmax= nzlmax;
+        lsub = m_Glu->lsub;
+      }
+      int ifrom, ito = nextl; 
+      for (ifrom = xlsub(jcol); ifrom < nextl;)
+        lsub(ito++) = lsub(ifrom++);
+      for (i = jcol+1; i <=kcol; i++)xlsub(i) = nextl;
+      nextl = ito;
+    }
+    xsup(nsuper+1) = kcol + 1; // Start of next available supernode
+    supno(kcol+1) = nsuper;
+    xprune(kcol) = nextl;
+    xlsub(kcol+1) = nextl;
+    return 0;
+  }
+   
+  
+#endif
\ No newline at end of file