From 4126cb636907246c4dab1987ace7eba0cb66473b Mon Sep 17 00:00:00 2001
From: Gael Guennebaud <g.gael@free.fr>
Date: Fri, 17 Oct 2014 16:52:56 +0200
Subject: [PATCH] Fix SparseLU::absDeterminant and add respective unit test
 (grafted from a370b1f2e20804f0b99bcda1c3940f00300e4f90 )

---
 Eigen/src/SparseLU/SparseLU.h                 |  5 ++-
 .../src/SparseLU/SparseLU_SupernodalMatrix.h  |  4 +--
 test/sparse_solver.h                          | 33 +++++++++++++++++++
 test/sparselu.cpp                             |  3 ++
 4 files changed, 40 insertions(+), 5 deletions(-)

diff --git a/Eigen/src/SparseLU/SparseLU.h b/Eigen/src/SparseLU/SparseLU.h
index 1d592f2c8..d6f5d59ca 100644
--- a/Eigen/src/SparseLU/SparseLU.h
+++ b/Eigen/src/SparseLU/SparseLU.h
@@ -260,14 +260,13 @@ class SparseLU : public internal::SparseLUImpl<typename _MatrixType::Scalar, typ
       eigen_assert(m_factorizationIsOk && "The matrix should be factorized first.");
       // Initialize with the determinant of the row matrix
       Scalar det = Scalar(1.);
-      //Note that the diagonal blocks of U are stored in supernodes,
+      // Note that the diagonal blocks of U are stored in supernodes,
       // which are available in the  L part :)
       for (Index j = 0; j < this->cols(); ++j)
       {
         for (typename SCMatrix::InnerIterator it(m_Lstore, j); it; ++it)
         {
-          if(it.row() < j) continue;
-          if(it.row() == j)
+          if(it.index() == j)
           {
             det *= (std::abs)(it.value());
             break;
diff --git a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
index ad6f2183f..b16afd6a4 100644
--- a/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
+++ b/Eigen/src/SparseLU/SparseLU_SupernodalMatrix.h
@@ -189,8 +189,8 @@ class MappedSuperNodalMatrix<Scalar,Index>::InnerIterator
         m_idval(mat.colIndexPtr()[outer]),
         m_startidval(m_idval),
         m_endidval(mat.colIndexPtr()[outer+1]),
-        m_idrow(mat.rowIndexPtr()[outer]),
-        m_endidrow(mat.rowIndexPtr()[outer+1])
+        m_idrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]]),
+        m_endidrow(mat.rowIndexPtr()[mat.supToCol()[mat.colToSup()[outer]]+1])
     {}
     inline InnerIterator& operator++()
     { 
diff --git a/test/sparse_solver.h b/test/sparse_solver.h
index d84aff070..244e81c1b 100644
--- a/test/sparse_solver.h
+++ b/test/sparse_solver.h
@@ -124,7 +124,23 @@ void check_sparse_determinant(Solver& solver, const typename Solver::MatrixType&
   Scalar refDet = dA.determinant();
   VERIFY_IS_APPROX(refDet,solver.determinant());
 }
+template<typename Solver, typename DenseMat>
+void check_sparse_abs_determinant(Solver& solver, const typename Solver::MatrixType& A, const DenseMat& dA)
+{
+  using std::abs;
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  
+  solver.compute(A);
+  if (solver.info() != Success)
+  {
+    std::cerr << "sparse solver testing: factorization failed (check_sparse_abs_determinant)\n";
+    return;
+  }
 
+  Scalar refDet = abs(dA.determinant());
+  VERIFY_IS_APPROX(refDet,solver.absDeterminant());
+}
 
 template<typename Solver, typename DenseMat>
 int generate_sparse_spd_problem(Solver& , typename Solver::MatrixType& A, typename Solver::MatrixType& halfA, DenseMat& dA, int maxSize = 300)
@@ -324,3 +340,20 @@ template<typename Solver> void check_sparse_square_determinant(Solver& solver)
     check_sparse_determinant(solver, A, dA);
   }
 }
+
+template<typename Solver> void check_sparse_square_abs_determinant(Solver& solver)
+{
+  typedef typename Solver::MatrixType Mat;
+  typedef typename Mat::Scalar Scalar;
+  typedef Matrix<Scalar,Dynamic,Dynamic> DenseMatrix;
+
+  // generate the problem
+  Mat A;
+  DenseMatrix dA;
+  generate_sparse_square_problem(solver, A, dA, 30);
+  A.makeCompressed();
+  for (int i = 0; i < g_repeat; i++) {
+    check_sparse_abs_determinant(solver, A, dA);
+  }
+}
+
diff --git a/test/sparselu.cpp b/test/sparselu.cpp
index 37980defc..52371cb12 100644
--- a/test/sparselu.cpp
+++ b/test/sparselu.cpp
@@ -44,6 +44,9 @@ 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_abs_determinant(sparselu_colamd);
+  check_sparse_square_abs_determinant(sparselu_amd);
 }
 
 void test_sparselu()