GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-08-12 19:59:05 +08:00
Code Issues Packages Projects Releases Wiki Activity

bug #792: SparseLU::factorize failed for structurally rank deficient matrices

Browse Source
This commit is contained in:
Christoph Hertzberg 2015-07-26 20:30:30 +02:00
parent 4b3052c54d
commit a44d022caf
3 changed files with 27 additions and 8 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
Eigen/src/SparseLU/SparseLU_pivotL.h
View File

@ -71,7 +71,7 @@ Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScal
// Determine the largest abs numerical value for partial pivoting // Determine the largest abs numerical value for partial pivoting
Index diagind = iperm_c(jcol); // diagonal index Index diagind = iperm_c(jcol); // diagonal index
RealScalar pivmax = 0.0; RealScalar pivmax(-1.0);
Index pivptr = nsupc; Index pivptr = nsupc;
Index diag = emptyIdxLU; Index diag = emptyIdxLU;
RealScalar rtemp; RealScalar rtemp;
@ -87,8 +87,9 @@ Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScal
} }
// Test for singularity // Test for singularity
if ( pivmax == 0.0 ) { if ( pivmax <= RealScalar(0.0) ) {
pivrow = lsub_ptr[pivptr]; // 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); perm_r(pivrow) = StorageIndex(jcol);
return (jcol+1); return (jcol+1);
} }
@ -104,13 +105,13 @@ Index SparseLUImpl<Scalar,StorageIndex>::pivotL(const Index jcol, const RealScal
// Diagonal element exists // Diagonal element exists
using std::abs; using std::abs;
rtemp = abs(lu_col_ptr[diag]); 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]; pivrow = lsub_ptr[pivptr];
} }
// Record pivot row // Record pivot row
perm_r(pivrow) = StorageIndex(jcol); perm_r(pivrow) = StorageIndex(jcol);
// Interchange row subscripts // Interchange row subscripts
if (pivptr != nsupc ) if (pivptr != nsupc )
{ {

20
test/sparse_solver.h
View File

@ -332,7 +332,18 @@ Index generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, De
return size; return size;
} }
template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000)
struct prune_column {
Index m_col;
prune_column(Index col) : m_col(col) {}
template<class Scalar>
bool operator()(Index, Index col, const Scalar&) const {
return col != m_col;
}
};
template<typename Solver> void check_sparse_square_solving(Solver& solver, int maxSize = 300, int maxRealWorldSize = 100000, bool checkDeficient = false)
{ {
typedef typename Solver::MatrixType Mat; typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar; typedef typename Mat::Scalar Scalar;
@ -364,6 +375,13 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver, int m
b = DenseVector::Zero(size); b = DenseVector::Zero(size);
check_sparse_solving(solver, A, b, dA, b); check_sparse_solving(solver, A, b, dA, b);
} }
// regression test for Bug 792 (structurally rank deficient matrices):
if(checkDeficient && size>1) {
Index col = internal::random<int>(0,size-1);
A.prune(prune_column(col));
solver.compute(A);
VERIFY_IS_EQUAL(solver.info(), NumericalIssue);
}
} }
// First, get the folder // First, get the folder

4
test/sparselu.cpp
View File

@ -42,8 +42,8 @@ template<typename T> void test_sparselu_T()
SparseLU<SparseMatrix<T, ColMajor, long int>, NaturalOrdering<long int> > sparselu_natural; SparseLU<SparseMatrix<T, ColMajor, long int>, NaturalOrdering<long int> > sparselu_natural;
check_sparse_square_solving(sparselu_colamd); check_sparse_square_solving(sparselu_colamd);
check_sparse_square_solving(sparselu_amd, 300, 2000); check_sparse_square_solving(sparselu_amd, 300, 2000, !true); // FIXME AMD ordering fails for structurally deficient matrices!
check_sparse_square_solving(sparselu_natural, 300, 2000); check_sparse_square_solving(sparselu_natural, 300, 2000, true);
check_sparse_square_abs_determinant(sparselu_colamd); check_sparse_square_abs_determinant(sparselu_colamd);
check_sparse_square_abs_determinant(sparselu_amd); check_sparse_square_abs_determinant(sparselu_amd);