bug #792: SparseLU::factorize failed for structurally rank deficient matrices

2 files changed, 21 insertions, 3 deletions

diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index 64f6ac4fe..90362ca0f 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -332,7 +332,18 @@ Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, De
   return size;
 }
 
-template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000)
+
+struct prune_column {
+  Index m_col;
+  prune_column(Index col) : m_col(col) {}
+  template<class Scalar>
+  bool operator()(Index, Index col, const Scalar&) const {
+    return col != m_col;
+  }
+};
+
+
+template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000, bool checkDeficient = false)
 {
   typedef typename Solver::MatrixType Mat;
   typedef typename Mat::Scalar Scalar;
@@ -364,6 +375,13 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver, int m
       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) {
+      Index 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 231c857ad..c725847d8 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -42,8 +42,8 @@ template<typename T> void test_sparselu_T()
   SparseLU<SparseMatrix<T, ColMajor, long int>, NaturalOrdering<long int> > sparselu_natural;
   
   check_sparse_square_solving(sparselu_colamd); 
-  check_sparse_square_solving(sparselu_amd, 300, 2000);
-  check_sparse_square_solving(sparselu_natural, 300, 2000);
+  check_sparse_square_solving(sparselu_amd, 300, 2000, !true); // FIXME AMD ordering fails for structurally deficient matrices!
+  check_sparse_square_solving(sparselu_natural, 300, 2000, true);
   
   check_sparse_square_abs_determinant(sparselu_colamd);
   check_sparse_square_abs_determinant(sparselu_amd);