aboutsummaryrefslogtreecommitdiffhomepage
path: root/Eigen/src/SparseLU
diff options
context:
space:
mode:
authorGravatar Christoph Hertzberg <chtz@informatik.uni-bremen.de>2015-07-26 20:30:30 +0200
committerGravatar Christoph Hertzberg <chtz@informatik.uni-bremen.de>2015-07-26 20:30:30 +0200
commita44d022cafce5b1f8e452a22d43bd23265438ca6 (patch)
tree7ebe6a90bdf3e4aa70554e8c8beccc3443e1fbb7 /Eigen/src/SparseLU
parent4b3052c54d14f486c2665530f06cdeff636169e4 (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.h11
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 )
{