Created attachment 646 [details]
The attached testcase crashes when compiled without NDEBUG on this assertion:
IncompleteCholesky.h:155: void Eigen::IncompleteCholesky<Scalar, _UpLo, _OrderingType>::_solve_impl(const Rhs&, Dest&) const [with Rhs = Eigen::Matrix<double, -1, 1>; Dest = Eigen::Matrix<double, -1, 1>; Scalar = double; int _UpLo = 1; _OrderingType = Eigen::AMDOrdering<int>]: Assertion `m_factorizationIsOk && "factorize() should be called first"' failed.
With N < 10, the assertion does not trigger.
This is because the solver failed, as reported by:
and more precisely, the preconditioner failed to factorize the matrix:
Of course, ILLT should not fail on a SPD matrix, so let's see what's going wrong.
OK, so we probably need to restart with a larger shift until the factorization succeed, as in algorithm 3.1 of: http://ftp.mcs.anl.gov/pub/tech_reports/reports/P682.pdf
Should be straightforward. Patch welcome ;)
Summary: Bug 1150: make IncompleteCholesky more robust by iteratively increase the shift until the factorization succeed (with at most 10 attempts).
(In reply to Gael Guennebaud from comment #3)
> Summary: Bug 1150: make IncompleteCholesky more robust by iteratively
> increase the shift until the factorization succeed (with at most 10
Thanks for fixing it. Is there still a risk of failure with the shift increment strategy ? If yes, would it occur on specific cases like big matricies, ill conditioned matricies, ... ?
Yes there is still a risk (e.g., with a large negative eigenvalue). It performs a maximum of 10 iteration, doubling the shift each time. If it still fails, you can increase the initial shift with the respective method, but with such a large shift, better use another preconditioning technique.