diff options
author | Christoph Hertzberg <chtz@informatik.uni-bremen.de> | 2015-07-26 20:30:30 +0200 |
---|---|---|
committer | Christoph Hertzberg <chtz@informatik.uni-bremen.de> | 2015-07-26 20:30:30 +0200 |
commit | a44d022cafce5b1f8e452a22d43bd23265438ca6 (patch) | |
tree | 7ebe6a90bdf3e4aa70554e8c8beccc3443e1fbb7 /Eigen/src/SparseLU | |
parent | 4b3052c54d14f486c2665530f06cdeff636169e4 (diff) |
bug #792: SparseLU::factorize failed for structurally rank deficient matrices
Diffstat (limited to 'Eigen/src/SparseLU')
-rw-r--r-- | Eigen/src/SparseLU/SparseLU_pivotL.h | 11 |
1 files changed, 6 insertions, 5 deletions
diff --git a/Eigen/src/SparseLU/SparseLU_pivotL.h b/Eigen/src/SparseLU/SparseLU_pivotL.h index 562128b69..a86dac93f 100644 --- a/Eigen/src/SparseLU/SparseLU_pivotL.h +++ b/Eigen/src/SparseLU/SparseLU_pivotL.h @@ -71,7 +71,7 @@ Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScal // 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,StorageIndex>::pivotL(const Index jcol, const RealScal } // 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) = StorageIndex(jcol); return (jcol+1); } @@ -104,13 +105,13 @@ Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScal // Diagonal element exists using std::abs; rtemp = abs(lu_col_ptr[diag]); - if (rtemp != 0.0 && rtemp >= thresh) pivptr = diag; + if (rtemp != RealScalar(0.0) && rtemp >= thresh) pivptr = diag; } pivrow = lsub_ptr[pivptr]; } // Record pivot row - perm_r(pivrow) = StorageIndex(jcol); + perm_r(pivrow) = StorageIndex(jcol); // Interchange row subscripts if (pivptr != nsupc ) { |