GitHub-Proxy/eigen
1
0
Fork 0
mirror of https://gitlab.com/libeigen/eigen.git synced 2025-07-09 14:41:49 +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:39:32 +02:00
parent c4432aad15
commit 6a3797f46f
3 changed files with 23 additions and 5 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -71,7 +71,7 @@ Index SparseLUImpl<Scalar,Index>::pivotL(const Index jcol, const RealScalar& dia
// 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,Index>::pivotL(const Index jcol, const RealScalar& dia
} }
// 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) = jcol; perm_r(pivrow) = jcol;
return (jcol+1); return (jcol+1);
} }

19
test/sparse_solver.h
View File

@ -277,7 +277,17 @@ int generate_sparse_square_problem(Solver&, typename Solver::MatrixType& A, Dens
return size; return size;
} }
template<typename Solver> void check_sparse_square_solving(Solver& solver)
struct prune_column {
int m_col;
prune_column(int col) : m_col(col) {}
template<class Scalar>
bool operator()(int, int col, const Scalar&) const {
return col != m_col;
}
};
template<typename Solver> void check_sparse_square_solving(Solver& solver, bool checkDeficient = false)
{ {
typedef typename Solver::MatrixType Mat; typedef typename Solver::MatrixType Mat;
typedef typename Mat::Scalar Scalar; typedef typename Mat::Scalar Scalar;
@ -309,6 +319,13 @@ template<typename Solver> void check_sparse_square_solving(Solver& solver)
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) {
int 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

2
test/sparselu.cpp
View File

@ -43,7 +43,7 @@ template<typename T> void test_sparselu_T()
check_sparse_square_solving(sparselu_colamd); check_sparse_square_solving(sparselu_colamd);
check_sparse_square_solving(sparselu_amd); check_sparse_square_solving(sparselu_amd);
check_sparse_square_solving(sparselu_natural); check_sparse_square_solving(sparselu_natural, 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);