Update Ordering interface

Eigen/src/OrderingMethods/Ordering.h

@@ -60,7 +60,9 @@ class AMDOrdering
   public:
     typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
     
-    /** Compute the permutation vector from a column-major sparse matrix */
+    /** Compute the permutation vector from a sparse matrix
+     * This routine is much faster if the input matrix is column-major     
+     */
     template <typename MatrixType>
     void operator()(const MatrixType& mat, PermutationType& perm)
     {
@@ -73,7 +75,7 @@ class AMDOrdering
       internal::minimum_degree_ordering(symm, perm);
     }
     
-    /** Compute the permutation with a self adjoint matrix */
+    /** Compute the permutation with a selfadjoint matrix */
     template <typename SrcType, unsigned int SrcUpLo> 
     void operator()(const SparseSelfAdjointView<SrcType, SrcUpLo>& mat, PermutationType& perm)
     { 
@@ -85,6 +87,26 @@ class AMDOrdering
     }
 };
 
+/** 
+ * Get the natural ordering
+ * 
+ *NOTE Returns an empty permutation matrix
+ * \tparam  Index The type of indices of the matrix 
+ */
+template <typename Index>
+class NaturalOrdering
+{
+  public:
+    typedef PermutationMatrix<Dynamic, Dynamic, Index> PermutationType;
+    
+    /** Compute the permutation vector from a column-major sparse matrix */
+    template <typename MatrixType>
+    void operator()(const MatrixType& mat, PermutationType& perm)
+    {
+      perm.resize(0); 
+    }
+    
+};
 
 /** 
  * Get the column approximate minimum degree ordering 

Eigen/src/SparseLU/SparseLU.h

@@ -255,7 +255,7 @@ class SparseLU
     void initperfvalues()
     {
       m_panel_size = 12; 
-      m_relax = 1; 
+      m_relax = 6; 
       m_maxsuper = 100; 
       m_rowblk = 200; 
       m_colblk = 60; 
@@ -320,26 +320,31 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
   // Compute the fill-reducing ordering  
   // TODO Currently, the only  available ordering method is AMD. 
   
-  OrderingType ord(mat); 
-  m_perm_c = ord.get_perm(); 
+  OrderingType ord; 
+  ord(mat,m_perm_c);
   //FIXME Check the right semantic behind m_perm_c
   // that is, column j of mat goes to column m_perm_c(j) of mat * m_perm_c; 
   
+  //DEBUG : Set the natural ordering 
+  for (int i = 0; i < mat.cols(); i++)
+    m_perm_c.indices()(i) = i; 
  
   // Apply the permutation to the column of the input  matrix
-  m_mat = mat * m_perm_c; 
+  m_mat = mat * m_perm_c.inverse(); 
   
     
   // Compute the column elimination tree of the permuted matrix 
   if (m_etree.size() == 0)  m_etree.resize(m_mat.cols());
+  
   LU_sp_coletree(m_mat, m_etree); 
-    
+     
   // In symmetric mode, do not do postorder here
   if (!m_symmetricmode) {
     IndexVector post, iwork; 
     // Post order etree
     LU_TreePostorder(m_mat.cols(), m_etree, post); 
       
+   
     // Renumber etree in postorder 
     int m = m_mat.cols(); 
     iwork.resize(m+1);
@@ -348,12 +353,15 @@ void SparseLU<MatrixType, OrderingType>::analyzePattern(const MatrixType& mat)
     
     // Postmultiply A*Pc by post, i.e reorder the matrix according to the postorder of the etree
     
-    PermutationType post_perm(m);; 
+    PermutationType post_perm(m); //FIXME Use vector constructor
     for (int i = 0; i < m; i++) 
       post_perm.indices()(i) = post(i); 
-    //m_mat = m_mat * post_perm; // FIXME This should surely be in factorize()  
+    
+//     m_mat = m_mat * post_perm.inverse(); // FIXME This should surely be in factorize()  
+    
     // Composition of the two permutations
     m_perm_c = m_perm_c * post_perm;
+    
   } // end postordering 
   
   m_analysisIsOk = true; 
@@ -402,9 +410,14 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   
   
   // Apply the column permutation computed in analyzepattern()
-  m_mat = matrix * m_perm_c; 
+  m_mat = matrix * m_perm_c.inverse(); 
   m_mat.makeCompressed(); 
   
+  // DEBUG ... Watch matrix permutation  
+  const int *asub_in = matrix.innerIndexPtr(); 
+  const int *colptr_in = matrix.outerIndexPtr(); 
+  int * asub = m_mat.innerIndexPtr(); 
+  int * colptr = m_mat.outerIndexPtr(); 
   int m = m_mat.rows();
   int n = m_mat.cols();
   int nnz = m_mat.nonZeros();
@@ -455,7 +468,8 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   
   // Setup Permutation vectors
   // Compute the inverse of perm_c
-  PermutationType iperm_c (m_perm_c.inverse() );
+//   PermutationType iperm_c (m_perm_c.inverse() );
+  PermutationType iperm_c (m_perm_c);
   
   // Identify initial relaxed snodes
   IndexVector relax_end(n);
@@ -464,6 +478,9 @@ void SparseLU<MatrixType, OrderingType>::factorize(const MatrixType& matrix)
   else
     LU_relax_snode<IndexVector>(n, m_etree, m_relax, marker, relax_end);
   
+    //DEBUG
+//   std::cout<< "relax_end " <<relax_end.transpose() << std::endl; 
+  
   m_perm_r.resize(m); 
   m_perm_r.indices().setConstant(-1); //FIXME
   marker.setConstant(-1);

Eigen/src/SparseLU/SparseLU_Coletree.h

@@ -75,7 +75,6 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
   IndexVector pp(nc); // disjoint sets 
   pp.setZero(); // Initialize disjoint sets 
   IndexVector firstcol(nr); // First nonzero column in each row 
-  firstcol.setZero(); 
   
   //Compute first nonzero column in each row 
   int row,col; 
@@ -95,7 +94,9 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
   int rset, cset, rroot; 
   for (col = 0; col < nc; col++) 
   {
-    cset = pp(col) = col; // Initially, each element is in its own set //FIXME
+//     cset = pp(col) = col; // Initially, each element is in its own set //FIXME
+    pp(col) = col; 
+    cset = col; 
     root(cset) = col; 
     parent(col) = nc; 
     for (typename MatrixType::InnerIterator it(mat, col); it; ++it)
@@ -107,7 +108,9 @@ int LU_sp_coletree(const MatrixType& mat, IndexVector& parent)
       if (rroot != col) 
       {
         parent(rroot) = col; 
-        cset = pp(cset) = rset; // Get the union of cset and rset  //FIXME
+//         cset = pp(cset) = rset; // Get the union of cset and rset  //FIXME
+        pp(cset) = rset; 
+        cset = rset; 
         root(cset) = col; 
       }
     }

Eigen/src/SparseLU/SparseLU_relax_snode.h

@@ -57,7 +57,7 @@ void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, Inde
 {
   
   // compute the number of descendants of each node in the etree
-  register int j, parent; 
+  int j, parent; 
   relax_end.setConstant(IND_EMPTY);
   descendants.setZero();
   for (j = 0; j < n; j++) 
@@ -66,9 +66,8 @@ void LU_relax_snode (const int n, IndexVector& et, const int relax_columns, Inde
     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 
+  int snode_start; // beginning of a snode 
   for (j = 0; j < n; )
   {
     parent = et(j);

Eigen/src/SparseLU/SparseLU_snode_dfs.h

@@ -68,8 +68,8 @@
     Index& nzlmax = glu.nzlmax; 
     int mem; 
     Index nsuper = ++supno(jcol); // Next available supernode number
-    register int nextl = xlsub(jcol); //Index of the starting location of the jcol-th column in lsub
-    register int i,k; 
+    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++)
     {
@@ -86,7 +86,7 @@
           if( nextl >= nzlmax )
           {
             mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu.num_expansions);
-            if (mem) return mem; 
+            if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
           }
         }
       }
@@ -100,7 +100,7 @@
       while (new_next > nzlmax)
       {
         mem = LUMemXpand<IndexVector>(lsub, nzlmax, nextl, LSUB, glu.num_expansions);
-        if (mem) return mem; 
+        if (mem) return mem; // Memory expansion failed... Return the memory allocated so far
       }
       Index ifrom, ito = nextl; 
       for (ifrom = xlsub(jcol); ifrom < nextl;)

Eigen/src/SuperLUSupport/SuperLUSupport.h

@@ -627,6 +627,9 @@ void SuperLU<MatrixType>::factorize(const MatrixType& a)
   
   this->initFactorization(a);
   
+  //DEBUG
+  m_sluOptions.ColPerm = NATURAL;
+  m_sluOptions.Equil = NO; 
   int info = 0;
   RealScalar recip_pivot_growth, rcond;
   RealScalar ferr, berr;

bench/spbench/test_sparseLU.cpp

@@ -17,7 +17,7 @@ int main(int argc, char **args)
   typedef Matrix<double, Dynamic, Dynamic> DenseMatrix;
   typedef Matrix<double, Dynamic, 1> DenseRhs;
   VectorXd b, x, tmp;
-  SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<double, int> >   solver;
+  SparseLU<SparseMatrix<double, ColMajor>, AMDOrdering<int> >   solver;
   ifstream matrix_file; 
   string line;
   int  n;
@@ -52,7 +52,7 @@ int main(int argc, char **args)
   }
 
   /* Compute the factorization */
-  solver.isSymmetric(true);
+//   solver.isSymmetric(true);
   solver.compute(A);
   
   solver._solve(b, x);