Improve the permutation

3 files changed, 39 insertions, 12 deletions

diff --git a/Eigen/src/SparseCore/SparseMatrix.h b/Eigen/src/SparseCore/SparseMatrix.h
index 214f130f5..52a9dab70 100644
--- a/Eigen/src/SparseCore/SparseMatrix.h
+++ b/Eigen/src/SparseCore/SparseMatrix.h
@@ -477,6 +477,18 @@ class SparseMatrix
       m_data.squeeze();
     }
 
+    /** Turns the matrix into the uncompressed mode */
+    void Uncompress()
+    {
+      if(m_innerNonZeros != 0)
+        return; 
+      m_innerNonZeros = new Index[m_outerSize]; 
+      for (int i = 0; i < m_outerSize; i++)
+      {
+        m_innerNonZeros[i] = m_outerIndex[i+1] - m_outerIndex[i]; 
+      }
+    }
+    
     /** Suppresses all nonzeros which are \b much \b smaller \b than \a reference under the tolerence \a epsilon */
     void prune(Scalar reference, RealScalar epsilon = NumTraits<RealScalar>::dummy_precision())
     {
diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 474dfdedc..70898958b 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -346,8 +346,17 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
 
   
   // Apply the permutation to the column of the input  matrix
-  m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix
-  
+//   m_mat = mat * m_perm_c.inverse(); //FIXME It should be less expensive here to permute only the structural pattern of the matrix
+   
+  //First copy the whole input matrix. 
+  m_mat = mat;
+  m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers. FIXME : This vector is filled but not subsequently used.  
+  //Then, permute only the column pointers
+  for (int i = 0; i < mat.cols(); i++)
+  {
+    m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i]; 
+    m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = mat.outerIndexPtr()[i+1] - mat.outerIndexPtr()[i]; 
+  }
     
   // Compute the column elimination tree of the permuted matrix 
   if (m_etree.size() == 0)  m_etree.resize(m_mat.cols());
@@ -424,8 +433,15 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   
   
   // Apply the column permutation computed in analyzepattern()
-  m_mat = matrix * m_perm_c.inverse(); 
-  m_mat.makeCompressed(); 
+  //   m_mat = matrix * m_perm_c.inverse(); 
+  m_mat = matrix;
+  m_mat.Uncompress(); //NOTE: The effect of this command is only to create the InnerNonzeros pointers.
+  //Then, permute only the column pointers
+  for (int i = 0; i < matrix.cols(); i++)
+  {
+    m_mat.outerIndexPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i]; 
+    m_mat.innerNonZeroPtr()[m_perm_c.indices()(i)] = matrix.outerIndexPtr()[i+1] - matrix.outerIndexPtr()[i]; 
+  }
   
   int m = m_mat.rows();
   int n = m_mat.cols();
@@ -504,7 +520,7 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
       
       // Factorize the relaxed supernode(jcol:kcol)
       // First, determine the union of the row structure of the snode 
-      info = LU_snode_dfs(jcol, kcol, m_mat.innerIndexPtr(), m_mat.outerIndexPtr(), xprune, marker, m_glu); 
+      info = LU_snode_dfs(jcol, kcol, m_mat, xprune, marker, m_glu); 
       if ( info ) 
       {
         std::cerr << "MEMORY ALLOCATION FAILED IN SNODE_DFS() \n";
diff --git a/Eigen/src/SparseLU/SparseLU_snode_dfs.h b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
index 6b2817262..150d9d0ef 100644
--- a/Eigen/src/SparseLU/SparseLU_snode_dfs.h
+++ b/Eigen/src/SparseLU/SparseLU_snode_dfs.h
@@ -57,8 +57,8 @@
  * \param marker (in/out) working vector
  * \return 0 on success, > 0 size of the memory when memory allocation failed
  */
-  template <typename IndexVector, typename ScalarVector>
-  int LU_snode_dfs(const int jcol, const int kcol, const typename IndexVector::Scalar* asub, const typename IndexVector::Scalar* colptr, IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
+  template <typename MatrixType, typename IndexVector, typename ScalarVector>
+  int LU_snode_dfs(const int jcol, const int kcol,const MatrixType& mat,  IndexVector& xprune, IndexVector& marker, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
   {
     typedef typename IndexVector::Scalar Index; 
     IndexVector& xsup = glu.xsup; 
@@ -69,14 +69,13 @@
     int mem; 
     Index nsuper = ++supno(jcol); // Next available supernode number
     int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
-    int i,k; 
     int krow,kmark; 
-    for (i = jcol; i <=kcol; i++)
+    for (int i = jcol; i <=kcol; i++)
     {
       // For each nonzero in A(*,i)
-      for (k = colptr[i]; k < colptr[i+1]; k++)
+      for (typename MatrixType::InnerIterator it(mat, i); it; ++it)
       {
-        krow = asub[k]; 
+        krow = it.row(); 
         kmark = marker(krow);
         if ( kmark != kcol )
         {
@@ -105,7 +104,7 @@
       Index ifrom, ito = nextl; 
       for (ifrom = xlsub(jcol); ifrom < nextl;)
         lsub(ito++) = lsub(ifrom++);
-      for (i = jcol+1; i <=kcol; i++) xlsub(i) = nextl;
+      for (int i = jcol+1; i <=kcol; i++) xlsub(i) = nextl;
       nextl = ito;
     }
     xsup(nsuper+1) = kcol + 1; // Start of next available supernode