Symbolic and numeric updates within the panel

Eigen/src/SparseLU/SparseLU.h

@@ -245,7 +245,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
   VectorXi& xlusup = m_GLu.xlusup;
   VectorXi& xusub = m_Glu.xusub;
     
-  supno(0) = -1; 
+  supno(0) = IND_EMPTY; 
   xsup(0) = xlsub(0) = xusub(0) = xlusup(0);
   int panel_size = m_panel_size; 
   int wdef = panel_size; // upper bound on panel width
@@ -262,7 +262,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
   int nseg; // Number of segments in each U-column 
   for (jcol = 0; jcol < min_mn; )
   {
-    if (relax_end(jcol) != -1) 
+    if (relax_end(jcol) != IND_EMPTY) 
     { // Starting a relaxed node from jcol
       kcol = relax_end(jcol); // End index of the relaxed snode 
       
@@ -298,7 +298,12 @@ void  SparseLU::factorize(const MatrixType& matrix)
         
         // Eliminate the current column 
         info = LU_pivotL(icol,  pivrow); 
-        eigen_assert(info == 0 && "The matrix is structurally singular"); 
+        if ( !info ) 
+        {
+          m_info = NumericalIssue; 
+          m_factorizationIsOk = false; 
+          return; 
+        }
       }
       jcol = icol; // The last column te be eliminated
     }
@@ -309,7 +314,7 @@ void  SparseLU::factorize(const MatrixType& matrix)
       panel_size = w_def;
       for (k = jcol + 1; k < std::min(jcol+panel_size, min_mn); k++)
       {
-        if (relax_end(k) != -1) 
+        if (relax_end(k) != IND_EMPTY) 
         {
           panel_size = k - jcol; 
           break; 
@@ -331,7 +336,9 @@ void  SparseLU::factorize(const MatrixType& matrix)
         
         nseg = nseg1; // begin after all the panel segments
         //Depth-first-search for the current column
-        info = LU_column_dfs(m, jj, ... ); 
+        VectorBlock<VectorXi> panel_lsubk(panel_lsub, k, m); //FIXME
+        VectorBlock<VectorXi> 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); 
         if ( !info ) 
         {
           m_info = NumericalIssue; 
@@ -339,7 +346,9 @@ void  SparseLU::factorize(const MatrixType& matrix)
           return; 
         }
         // Numeric updates to this column 
-        info = LU_column_bmod(jj, ... ); 
+        VectorBlock<VectorXi> dense_k(dense, k, m); //FIXME 
+        VectorBlock<VectorXi> 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); 
         if ( !info ) 
         {
           m_info = NumericalIssue; 
@@ -347,6 +356,20 @@ void  SparseLU::factorize(const MatrixType& matrix)
           return; 
         }
         
+        // Copy the U-segments to ucol(*)
+        
+        
+        // Form the L-segment 
+        info = LU_pivotL(...);
+        if ( !info ) 
+        {
+          m_info = NumericalIssue; 
+          m_factorizationIsOk = false; 
+          return; 
+        }
+        
+        // Prune columns (0:jj-1) using column jj
+        
       } // end for 
       jcol += panel_size;  // Move to the next panel
     } // end else 

Eigen/src/SparseLU/SparseLU_Memory.h

@@ -86,7 +86,7 @@ int SparseLU::LUMemInit(int lwork)
     nzlmax  = std::max(1, m_fill_ratio/4.) * annz; //???
   
     // Return the estimated size to the user if necessary
-    if (lwork = -1) 
+    if (lwork == IND_EMPTY) 
     {
       estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
                       + (nzlmax + nzumax) * iword + (nzlumax+nzumax) * dword + n); 
@@ -130,7 +130,7 @@ int SparseLU::LUMemInit(int lwork)
   }
   else  // m_fact == SamePattern_SameRowPerm;
   {
-    if (lwork = -1) 
+    if (lwork == IND_EMPTY) 
     {
       estimated_size = LU_GluIntArray(n) * iword + LU_TempSpace(m, m_panel_size)
                       + (Glu.nzlmax + Glu.nzumax) * iword + (Glu.nzlumax+Glu.nzumax) * dword + n); 
@@ -232,7 +232,7 @@ DestType* SparseLU::LUMemXpand(int jcol, int next, MemType mem_type, int& maxlen
     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");
+  eigen_assert(new_mem && "Can't expand memory"); // FIXME Should be an exception 
   
   return new_mem;
   

Eigen/src/SparseLU/SparseLU_Utils.h

@@ -28,5 +28,5 @@
 // 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)
+#define IND_EMPTY (-1)
 #endif

Eigen/src/SparseLU/SparseLU_column_bmod.h

@@ -0,0 +1,216 @@
+// 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 xcolumn_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_COLUMN_BMOD_H
+#define SPARSELU_COLUMN_BMOD_H
+/**
+ * \brief Performs numeric block updates (sup-col) in topological order
+ * 
+ * \param jcol current column to update
+ * \param nseg Number of segments in the U part
+ * \param dense Store the full representation of the column
+ * \param tempv working array 
+ * \param segrep segment representative ...
+ * \param repfnz ??? First nonzero column in each row ???  ...
+ * \param fpanelc First column in the current panel
+ * \param Glu Global LU data. 
+ * \return 0 - successful return 
+ *         > 0 - number of bytes allocated when run out of space
+ * 
+ */
+template <typename VectorType>
+int SparseLU::LU_column_bmod(const int jcol, const int nseg, VectorType& dense, VectorType& tempv, VectorXi& segrep, VectorXi& repfnz, int fpanelc, LU_GlobalLu_t& Glu)
+{
+    
+  int jsupno, k, ksub, krep, krep_ind, ksupno; 
+  /* krep = representative of current k-th supernode
+    * fsupc =  first supernodal column
+    * nsupc = number of columns in a supernode
+    * nsupr = number of rows in a supernode
+    * luptr = location of supernodal LU-block in storage
+    * kfnz = first nonz in the k-th supernodal segment
+    * no-zeros = no lf leading zeros in a supernodal U-segment
+    */
+  VectorXi& xsup = Glu.xsup; 
+  VectorXi& supno = Glu.supno; 
+  VectorXi& lsub = Glu.lsub; 
+  VectorXi& xlsub = Glu.xlsub; 
+  VectorXi& xlusup = Glu.xlusup; 
+  VectorType& lusup = Glu.lusup; 
+  int nzlumax = GLu.nzlumax; 
+  int jsupno = supno(jcol);
+  // For each nonzero supernode segment of U[*,j] in topological order 
+  k = nseg - 1; 
+  for (ksub = 0; ksub < nseg; ksub++)
+  {
+    krep = segrep(k); k--; 
+    ksupno = supno(krep); 
+    if (jsupno != ksupno )
+    {
+      // outside the rectangular supernode 
+      fsupc = xsup(ksupno); 
+      fst_col = std::max(fsupc, fpanelc); 
+      
+      // Distance from the current supernode to the current panel; 
+      // d_fsupc = 0 if fsupc > fpanelc
+      d_fsupc = fst_col - fsupc; 
+      
+      luptr = xlusup(fst_col) + d_fsupc; 
+      lptr = xlsub(fsupc) + d_fsupc; 
+      
+      kfnz = repfnz(krep); 
+      kfnz = std::max(kfnz, fpanelc); 
+      
+      segsize = krep - kfnz + 1; 
+      nsupc = krep - fst_col + 1; 
+      nsupr = xlsub(fsupc+1) - xlsub(fsupc); 
+      nrow = nsupr - d_fsupc - nsupc; 
+      krep_ind = lptr + nsupc - 1; 
+      
+      // NOTE  Unlike the original implementation in SuperLU, the only feature  
+      //  here is a sup-col update. 
+      
+      // Perform a triangular solver and block update, 
+      // then scatter the result of sup-col update to dense
+      no_zeros = kfnz - fst_col; 
+      // First, copy U[*,j] segment from dense(*) to tempv(*)
+      isub = lptr + no_zeros; 
+      for (i = 0; i ww segsize; i++)
+      {
+        irow = lsub(isub); 
+        tempv(i) = densee(irow); 
+        ++isub; 
+      }
+      // Dense triangular solve -- start effective triangle
+      luptr += nsupr * no_zeros + no_zeros; 
+      // Form Eigen matrix and vector 
+      Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) );
+      Map<VectorType> u(tempv.data(), segsize);
+      u = A.triangularView<Lower>().solve(u); 
+      
+      // Dense matrix-vector product y <-- A*x 
+      luptr += segsize; 
+      new (&A) (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); 
+      Map<VectorType> l( &(tempv.data()[segsize]), segsize); 
+      l= A * u;
+      
+      // Scatter tempv[] into SPA dense[] as a temporary storage 
+      isub = lptr + no_zeros; 
+      for (i = 0; i w segsize; i++)
+      {
+        irow = lsub(isub); 
+        dense(irow) = tempv(i); 
+        tempv(i) =  Scalar(0.0); 
+        ++isub;
+      }
+      
+      // Scatter l into SPA dense[]
+      for (i = 0; i < nrow; i++)
+      {
+        irow = lsub(isub); 
+        dense(irow) -= tempv(segsize + i); 
+        tempv(segsize + i) = Scalar(0.0); 
+        ++isub; 
+      }
+    } // end if jsupno 
+  } // end for each segment
+  
+  // Process the supernodal portion of  L\U[*,j]
+  nextlu = xlusup(jcol); 
+  fsupc = xsup(jsupno);
+  
+  // copy the SPA dense into L\U[*,j]
+  new_next = nextlu + xlsub(fsupc + 1) - xlsub(fsupc); 
+  while (new_next > nzlumax )
+  {
+    Glu.lusup = LUmemXpand<Scalar>(jcol, nextlu, LUSUP, &nzlumax); 
+    Glu.nzlumax = nzlumax; 
+    lusup = Glu.lusup; 
+    lsub = Glu.lsub; 
+  }
+  
+  for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
+  {
+    irow = lsub(isub);
+    lusub(nextlu) = dense(irow);
+    dense(irow) = Scalar(0.0); 
+    ++nextlu; 
+  }
+  
+  xlusup(jcol + 1) = nextlu;  // close L\U(*,jcol); 
+  
+  /* For more updates within the panel (also within the current supernode),
+   * should start from the first column of the panel, or the first column
+   * of the supernode, whichever is bigger. There are two cases:
+   *  1) fsupc < fpanelc, then fst_col <- fpanelc
+   *  2) fsupc >= fpanelc, then fst_col <-fsupc
+   */
+  fst_col = std::max(fsupc, fpanelc); 
+  
+  if (fst_col  < jcol)
+  {
+    // Distance between the current supernode and the current panel
+    // d_fsupc = 0 if fsupc >= fpanelc
+    d_fsupc = fst_col - fsupc; 
+    
+    lptr = xlsub(fsupc) + d_fsupc; 
+    luptr = xlusup(fst_col) + d_fsupc; 
+    nsupr = xlsub(fsupc+1) - xlsub(fsupc); // leading dimension
+    nsupc = jcol - fst_col; // excluding jcol 
+    nrow = nsupr - d_fsupc - nsupc; 
+    
+    // points to the beginning of jcol in snode L\U(jsupno) 
+    ufirst = xlusup(jcol) + d_fsupc; 
+    Map<Matrix<Scalar,Dynamic,Dynamic>, 0,  OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) ); 
+    Map<VectorType> l( &(lusup.data()[ufirst]), nsupc ); 
+    u = A.triangularView().solve(u); 
+    
+    new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); 
+    Map<VectorType> l( &(lusup.data()[ufirst+nsupc]), nsupr ); 
+    l = l - A * u;
+    
+  } // End if fst_col
+  return 0; 
+}
+#endif

Eigen/src/SparseLU/SparseLU_column_dfs.h

@@ -0,0 +1,269 @@
+// 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 xcolumn_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_COLUMN_DFS_H
+#define SPARSELU_COLUMN_DFS_H
+/**
+ * \brief Performs a symbolic factorization on column jcol and decide the supernode boundary
+ * 
+ * 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. 
+ * The location of the first nonzero in each supernodal segment 
+ * (supernodal entry location) is also returned. 
+ * 
+ * \param m number of rows in the matrix
+ * \param jcol Current column 
+ * \param perm_r Row permutation
+ * \param [in,out] nseg Number of segments in current U[*,j] - new segments appended
+ * \param lsub_col defines the rhs vector to start the dfs
+ * \param [in,out] segrep Segment representatives - new segments appended 
+ * \param repfnz
+ * \param xprune 
+ * \param marker
+ * \param parent
+ * \param xplore
+ * \param Glu global LU data 
+ * \return 0 success
+ *         > 0 number of bytes allocated when run out of space
+ * 
+ */
+int SparseLU::LU_column_dfs(const int m, const int jcol, VectorXi& perm_r, VectorXi& nseg  VectorXi& lsub_col, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu)
+{
+  typedef typename VectorXi::Index; 
+  
+  int jcolp1, jcolm1, jsuper, nsuper, nextl; 
+  int krow; // Row index of the current element 
+  int kperm; // permuted row index
+  int krep; // Supernode reprentative of the current row
+  int k, kmark; 
+  int chperm, chmark, chrep, oldrep, kchild; 
+  int myfnz; // First nonzero element in the current column
+  int xdfs, maxdfs, kpar;
+  
+  // Initialize pointers 
+  VectorXi& xsup = Glu.xsup; 
+  VectorXi& supno = Glu.supno; 
+  VectorXi& lsub = Glu.lsub; 
+  VectorXi& xlsub = Glu.xlsub; 
+  
+  nsuper = supno(jcol); 
+  jsuper = nsuper; 
+  nextl = xlsup(jcol); 
+  VectorBlock<VectorXi> marker2(marker, 2*m, m); 
+  // For each nonzero in A(*,jcol) do dfs 
+  for (k = 0; lsub_col[k] != IND_EMPTY; k++) 
+  {
+    krow = lsub_col(k); 
+    lsub_col(k) = IND_EMPTY; 
+    kmark = marker2(krow); 
+    
+    // krow was visited before, go to the next nonz; 
+    if (kmark == jcol) continue; 
+    
+    // For each unmarker nbr krow of jcol
+    //  krow is in L: place it in structure of L(*,jcol)
+    marker2(krow) = jcol; 
+    kperm = perm_r(krow); 
+    
+    if (kperm == IND_EMPTY ) 
+    {
+      lsub(nextl++) = krow; // krow is indexed into A
+      if ( nextl >= nzlmax )
+      {
+        Glu.lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax); 
+        //FIXME try... catch out of space 
+        Glu.nzlmax = nzlmax; 
+        lsub = Glu.lsub; 
+      }
+      if (kmark != jcolm1) jsuper = IND_EMPTY; // Row index subset testing
+    }
+    else 
+    {
+      // krow is in U : if its supernode-rep krep
+      // has been explored, update repfnz(*)
+      krep = xsup(supno(kperm)+1) - 1;
+      myfnz = repfnz(krep); 
+      
+      if (myfnz != IND_EMPTY )
+      {
+        // visited before 
+        if (myfnz > kperm) repfnz(krep) = kperm; 
+        // continue; 
+      }
+      else 
+      {
+        // otherwise, perform dfs starting at krep
+        oldrep = IND_EMPTY; 
+        parent(krep) = oldrep; 
+        repfnz(krep) = kperm; 
+        xdfs = xlsub(krep); 
+        maxdfs = xprune(krep); 
+        
+        do
+        {
+          // For each unmarked kchild of krep 
+          while (xdfs < maxdfs)
+          {
+            kchild = lsub(xdfs); 
+            xdfs++; 
+            chmark = marker2(kchild); 
+            
+            if (chmark != jcol) 
+            {
+             // Not reached yet 
+              marker2(kchild) = jcol; 
+              chperm = perm_r(kchild); 
+              
+              // if kchild is in L: place it in L(*,k)
+              if (chperm == IND_EMPTY)
+              {
+                lsub(nextl++) = kchild; 
+                if (nextl >= nzlmax)
+                {
+                  Glu.lsub = LUMemXpand<Index>(jcol, nextl, LSUB, nzlmax); 
+                  //FIXME Catch out of space errors 
+                  GLu.nzlmax = nzlmax; 
+                  lsub = Glu.lsub; 
+                }
+                if (chmark != jcolm1) jsuper = IND_EMPTY; 
+              } 
+              else 
+              {
+                // if 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(chrep); 
+                if (myfnz != IND_EMPTY) 
+                {
+                  // Visited before 
+                  if ( myfnz > chperm) repfnz(chrep) = chperm; 
+                }
+                else 
+                {
+                  // continue dfs at super-rep of kchild 
+                  xplore(krep) = xdfs; 
+                  oldrep = krep; 
+                  krep = chrep; // Go deeped down G(L^t)
+                  parent(krep) = olddrep; 
+                  repfnz(krep) = chperm; 
+                  xdfs = xlsub(krep); 
+                  maxdfs = xprune(krep); 
+                } // else myfnz 
+              } // else for chperm 
+              
+            } // if chmark 
+            
+          } // end while 
+          
+          // krow has no more unexplored nbrs; 
+          // place supernode-rep krep in postorder DFS.
+          // backtrack dfs to its parent
+          
+          segrep(nseg) = ;krep; 
+          ++nseg; 
+          kpar = parent(krep); // Pop from stack, mimic recursion
+          if (kpar == IND_EMPTY) break; // dfs done 
+          krep = kpar; 
+          xdfs = xplore(krep); 
+          maxdfs = xprune(krep); 
+          
+        } while ( kpar != IND_EMPTY); 
+        
+      } // else myfnz
+      
+    } // else kperm 
+    
+  } // for each nonzero ... 
+  
+  // check to see if j belongs in the same supeprnode as j-1
+  if ( jcol == 0 )
+  { // Do nothing for column 0 
+    nsuper = supno(0) = 0 ;
+  }
+  else 
+  {
+    fsupc = xsup(nsuper); 
+    jptr = xlsub(jcol); // Not yet compressed
+    jm1ptr = xlsub(jcolm1); 
+    
+    // Make sure the number of columns in a supernode doesn't
+    // exceed threshold
+    if ( (jcol - fsupc) >= m_maxsuper) jsuper = IND_EMPTY; 
+    
+    /* If jcol starts a new supernode, reclaim storage space in
+     * lsub from previous supernode. Note we only store 
+     * the subscript set of the first and last columns of 
+     * a supernode. (first for num values, last for pruning)
+     */
+    if (jsuper == IND_EMPTY)
+    { // starts a new supernode 
+      if ( (fsupc < jcolm1-1) ) 
+      { // >= 3 columns in nsuper
+        ito = xlsub(fsupcc+1)
+        xlsub(jcolm1) = ito; 
+        istop = ito + jptr - jm1ptr; 
+        xprune(jcolm1) = istop; // intialize xprune(jcol-1)
+        xlsub(jcol) = istop; 
+        
+        for (ifrom = jm1ptr; ifrom < nextl; ++ifrom, ++ito)
+          lsub(ito) = lsub(ifrom); 
+        nextl = ito;  // = istop + length(jcol)
+      }
+      nsuper++; 
+      supno(jcol) = nsuper; 
+    } // if a new supernode 
+  } // end else:  jcol > 0
+  
+  // Tidy up the pointers before exit
+  xsup(nsuper+1) = jcolp1; 
+  supno(jcolp1) = nsuper; 
+  xprune(jcol) = nextl;  // Intialize upper bound for pruning
+  xlsub(jcolp1) = nextl; 
+  
+  return 0; 
+}
+#endif

Eigen/src/SparseLU/SparseLU_panel_bmod.h

@@ -24,7 +24,7 @@
 
 /* 
  
- * NOTE: This file is the modified version of xpanel_dfs.c file in SuperLU 
+ * NOTE: This file is the modified version of xpanel_bmod.c file in SuperLU 
  
  * -- SuperLU routine (version 3.0) --
  * Univ. of California Berkeley, Xerox Palo Alto Research Center,
@@ -111,7 +111,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
       VectorBLock<VectorXi> dense_col(dense.segment(nextl_col, m)); // Scatter/gather entire matrix column from/to here
       
       kfnz = repfnz_col(krep); 
-      if ( kfnz == -1 ) 
+      if ( kfnz == IND_EMPTY ) 
         continue; // skip any zero segment
       
       segsize = krep - kfnz + 1;
@@ -143,7 +143,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
       
       luptr += segsize; 
       // Dense Matrix vector product y <-- A*x; 
-      new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), segsize, segsize, OuterStride<>(nsupr) ); 
+      new (&A) Map<Matrix<Scalar,Dynamic, Dynamic>, 0, OuterStride<> > ( &(lusup.data()[luptr]), nrow, segsize, OuterStride<>(nsupr) ); 
       Map<VectorType> l( &(tempv.data()[segsize]), segsize); 
       l= A * u;
       
@@ -157,7 +157,7 @@ void SparseLU::LU_panel_bmod(const int m, const int w, const int jcol, const int
       {
         irow = lsub(isub); 
         dense_col(irow) = tempv(i);
-        tempv(i) = zero; 
+        tempv(i) = Scalar(0.0); 
         isub++;
       }
      

Eigen/src/SparseLU/SparseLU_panel_dfs.h

@@ -72,7 +72,7 @@
  * 
  */
 template <typename MatrixType, typename VectorType>
-void 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)
+void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType& A, VectorXi& perm_r, int& nseg, VectorType& dense,  VectorXi& panel_lsub, VectorXi& segrep, VectorXi& repfnz, VectorXi& xprune, VectorXi& marker, VectorXi& parent, VectorXi& xplore, LU_GlobalLu_t& Glu)
 {
   
   int jj; // Index through each column in the panel 
@@ -115,7 +115,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
       // For each unmarked krow of jj
       marker(krow) = jj; 
       kperm = perm_r(krow); 
-      if (kperm == -1 ) {
+      if (kperm == IND_EMPTY ) {
         // krow is in L : place it in structure of L(*, jj)
         panel_lsub(nextl_col++) = krow;  // krow is indexed into A
       }
@@ -126,7 +126,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
         krep = xsup(supno(kperm)+1) - 1; 
         myfnz = repfnz_col(krep); 
         
-        if (myfnz != -1 )
+        if (myfnz != IND_EMPTY )
         {
           // Representative visited before
           if (myfnz > kperm ) repfnz_col(krep) = kperm; 
@@ -135,7 +135,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
         else 
         {
           // Otherwise, perform dfs starting at krep
-          oldrep = -1; 
+          oldrep = IND_EMPTY; 
           parent(krep) = oldrep; 
           repfnz_col(krep) = kperm; 
           xdfs =  xlsub(krep); 
@@ -155,7 +155,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
                 marker(kchild) = jj; 
                 chperm = perm_r(kchild); 
                 
-                if (chperm == -1) 
+                if (chperm == IND_EMPTY) 
                 {
                   // case kchild is in L: place it in L(*, j)
                   panel_lsub(nextl_col++) = kchild; 
@@ -168,7 +168,7 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
                   chrep = xsup(supno(chperm)+1) - 1; 
                   myfnz = repfnz_col(chrep); 
                   
-                  if (myfnz != -1) 
+                  if (myfnz != IND_EMPTY) 
                   { // Visited before 
                     if (myfnz > chperm) 
                       repfnz_col(chrep) = chperm; 
@@ -202,13 +202,13 @@ void SparseLU::LU_panel_dfs(const int m, const int w, const int jcol, MatrixType
             }
             
             kpar = parent(krep); // Pop recursion, mimic recursion 
-            if (kpar == -1) 
+            if (kpar == IND_EMPTY) 
               break; // dfs done 
             krep = kpar; 
             xdfs = xplore(krep); 
             maxdfs = xprune(krep); 
 
-          } while (kpar != -1); // Do until empty stack 
+          } while (kpar != IND_EMPTY); // Do until empty stack 
           
         } // end if (myfnz = -1)
 

Eigen/src/SparseLU/SparseLU_snode_bmod.h

@@ -54,9 +54,9 @@ int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc,
   VectorXi& 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; 
+  int irow, isub; 
   // Process the supernodal portion of L\U[*,jcol]
-  for (int isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
+  for (isub = xlsub(fsupc); isub < xlsub(fsupc+1); isub++)
   {
     irow = lsub(isub); 
     lusup(nextlu) = dense(irow);
@@ -72,20 +72,16 @@ int SparseLU::LU_dsnode_bmod (const int jcol, const int jsupno, const int fsupc,
     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.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) );
-    Map<Matrix<Scalar,Dynamic,Dynamic>, 0, OuterStride<> > A( &(lusup.data()[luptr]), nsupc, nsupc, OuterStride<>(nsupr) );
+    Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst]), nsupc);
-    Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst]), nsupc);
+    u = A.triangularView<Lower>().solve(u);
-    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)
-    Map<Matrix<Scalar,Dynamic,1> > u(&(lusup.data()[ufirst+nsupc], nsupc); 
+    new (&A) Map<Matrix<Scalar,Dynamic,Dynamic>,0,OuterStride<> > ( &(lusup.data()[luptr+nsupc]), nrow, nsupc, OuterStride<>(nsupr) ); 
-    u = A * l; 
+    Map<Matrix<Scalar,Dynamic,1> > l(&(lusup.data()[ufirst+nsupc], nsupc); 
-//     BLASFUNC(gemv)("N", &nrow, &nsupc, &alpha, &lusup[luptr+nsupc], &nsupr, &lusup[ufirst], &incx, &beta, &lusup[ufirst+nsupc], &incy);
+    l = l - A * u; 
     
     return 0;
 }