From 6a3797f46f6f02bb75c4d9df4c46796a733e14e2 Mon Sep 17 00:00:00 2001
From: Christoph Hertzberg <chtz@informatik.uni-bremen.de>
Date: Sun, 26 Jul 2015 20:39:32 +0200
Subject: [PATCH] bug #792: SparseLU::factorize failed for structurally rank
 deficient matrices

---
 Eigen/src/SparseLU/SparseLU_pivotL.h |  7 ++++---
 test/sparse_solver.h                 | 19 ++++++++++++++++++-
 test/sparselu.cpp                    |  2 +-
 3 files changed, 23 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h
index 457789c78..2e49ef667 100644
--- a/Eigen/src/SparseLU/SparseLU_pivotL.h
+++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
@@ -71,7 +71,7 @@ Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& dia
   
   // Determine the largest abs numerical value for partial pivoting 
   Index diagind = iperm_c(jcol); // diagonal index 
-  RealScalar pivmax = 0.0; 
+  RealScalar pivmax(-1.0);
   Index pivptr = nsupc; 
   Index diag = emptyIdxLU; 
   RealScalar rtemp;
@@ -87,8 +87,9 @@ Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& dia
   }
   
   // Test for singularity
-  if ( pivmax == 0.0 ) {
-    pivrow = lsub_ptr[pivptr];
+  if ( pivmax <= RealScalar(0.0) ) {
+    // if pivmax == -1, the column is structurally empty, otherwise it is only numerically zero
+    pivrow = pivmax < RealScalar(0.0) ? diagind : lsub_ptr[pivptr];
     perm_r(pivrow) = jcol;
     return (jcol+1);
   }
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 59d77daa2..833c0d889 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -277,7 +277,17 @@ int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, Dens
   return size;
 }
 
-template<typename Solver> void check_sparse_square_solving(Solver& solver)
+
+struct prune_column {
+  int m_col;
+  prune_column(int col) : m_col(col) {}
+  template<class Scalar>
+  bool operator()(int, int col, const Scalar&) const {
+    return col != m_col;
+  }
+};
+
+template<typename Solver> void check_sparse_square_solving(Solver& solver, bool checkDeficient = false)
 {
   typedef typename Solver::MatrixType Mat;
   typedef typename Mat::Scalar Scalar;
@@ -309,6 +319,13 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
       b = DenseVector::Zero(size);
       check_sparse_solving(solver, A, b, dA, b);
     }
+    // regression test for Bug 792 (structurally rank deficient matrices):
+    if(checkDeficient && size>1) {
+      int col = internal::random<int>(0,size-1);
+      A.prune(prune_column(col));
+      solver.compute(A);
+      VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
+    }
   }
   
   // First, get the folder 
diff --git a/test/sparselu.cpp b/test/sparselu.cpp
index 37eb069a9..f57108fa9 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -43,7 +43,7 @@ template<typename T> void test_sparselu_T()
   
   check_sparse_square_solving(sparselu_colamd); 
   check_sparse_square_solving(sparselu_amd);
-  check_sparse_square_solving(sparselu_natural);
+  check_sparse_square_solving(sparselu_natural, true);
   
   check_sparse_square_abs_determinant(sparselu_colamd);
   check_sparse_square_abs_determinant(sparselu_amd);