From 925ace196c182759026d3eb3edc06565ab5f01ee Mon Sep 17 00:00:00 2001
From: Desire NUENTSA <desire.nuentsa_wakam@inria.fr>
Date: Thu, 19 Jul 2012 18:15:23 +0200
Subject: [PATCH] correct bug in the complex version

---
 Eigen/src/SparseLU/SparseLU_pivotL.h | 11 ++++++-----
 bench/spbench/test_sparseLU.cpp      |  3 ++-
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h
index 39151f1e0..0c767c23a 100644
--- a/Eigen/src/SparseLU/SparseLU_pivotL.h
+++ b/Eigen/src/SparseLU/SparseLU_pivotL.h
@@ -71,7 +71,8 @@ template <typename IndexVector, typename ScalarVector>
 int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivotthresh, IndexVector& perm_r, IndexVector& iperm_c, int& pivrow, LU_GlobalLU_t<IndexVector, ScalarVector>& glu)
 {
   typedef typename IndexVector::Scalar Index; 
-  typedef typename ScalarVector::Scalar Scalar; 
+  typedef typename ScalarVector::Scalar Scalar;
+  typedef typename ScalarVector::RealScalar RealScalar;  
   // Initialize pointers 
   IndexVector& lsub = glu.lsub; // Compressed row subscripts of L rectangular supernodes.
   IndexVector& xlsub = glu.xlsub; // pointers to the beginning of each column subscript in lsub
@@ -88,10 +89,10 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
   
   // Determine the largest abs numerical value for partial pivoting 
   Index diagind = iperm_c(jcol); // diagonal index 
-  Scalar pivmax = 0.0; 
+  RealScalar pivmax = 0.0; 
   Index pivptr = nsupc; 
   Index diag = IND_EMPTY; 
-  Scalar rtemp;
+  RealScalar rtemp;
   Index isub, icol, itemp, k; 
   for (isub = nsupc; isub < nsupr; ++isub) {
     rtemp = std::abs(lu_col_ptr[isub]);
@@ -109,7 +110,7 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
     return (jcol+1);
   }
   
-  Scalar thresh = diagpivotthresh * pivmax; 
+  RealScalar thresh = diagpivotthresh * pivmax; 
   
   // Choose appropriate pivotal element 
   
@@ -119,7 +120,7 @@ int LU_pivotL(const int jcol, const typename ScalarVector::RealScalar diagpivott
     {
       // Diagonal element exists
       rtemp = std::abs(lu_col_ptr[diag]);
-      if (rtemp != Scalar(0.0) && rtemp >= thresh) pivptr = diag;
+      if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag;
     }
     pivrow = lsub_ptr[pivptr];
   }
diff --git a/bench/spbench/test_sparseLU.cpp b/bench/spbench/test_sparseLU.cpp
index 31273add5..08b6c926e 100644
--- a/bench/spbench/test_sparseLU.cpp
+++ b/bench/spbench/test_sparseLU.cpp
@@ -14,6 +14,7 @@ using namespace Eigen;
 int main(int argc, char **args)
 {
   typedef complex<double> scalar; 
+//   typedef double scalar; 
   SparseMatrix<scalar, ColMajor> A; 
   typedef SparseMatrix<scalar, ColMajor>::Index Index;
   typedef Matrix<scalar, Dynamic, Dynamic> DenseMatrix;
@@ -34,7 +35,7 @@ int main(int argc, char **args)
   bool iscomplex=false, isvector=false;
   int sym;
   getMarketHeader(args[1], sym, iscomplex, isvector);
-  if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
+//   if (iscomplex) { cout<< " Not for complex matrices \n"; return -1; }
   if (isvector) { cout << "The provided file is not a matrix file\n"; return -1;}
   if (sym != 0) { // symmetric matrices, only the lower part is stored
     SparseMatrix<scalar, ColMajor> temp;