diff options
author | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-07-19 18:15:23 +0200 |
---|---|---|
committer | Desire NUENTSA <desire.nuentsa_wakam@inria.fr> | 2012-07-19 18:15:23 +0200 |
commit | 925ace196c182759026d3eb3edc06565ab5f01ee (patch) | |
tree | 0095b8e8e36dfcdab19609b3564d045213cdd981 | |
parent | 59642da88bf83709e918667680e4ed63af4c31e5 (diff) |
correct bug in the complex version
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_pivotL.h | 11 | ||||
-rw-r--r-- | 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; |