Crash in IncompleteCholesky with zeros on diagonal

Summary

Hi! Using the IncompleteCholesky preconditioner on a large matrix, I get a crash on a range assertion. In release mode, this results in memory corruption.

Environment

• Operating System : Linux Ubuntu
• Architecture : x64 ...
• Eigen Version : 3.3.7 (package) ichol_crash.mat.gzand latest
• Compiler Version : Gcc9.4
• Compile Flags : any
• Vector Extension : avx2 available

Minimal Example

Not really a minimal example but a matrix market save is attached.

#include <eigen3/Eigen/Sparse>
#include <eigen3/unsupported/Eigen/SparseExtra>

int main() {
    Eigen::SparseMatrix<float> mat;
    Eigen::loadMarket(mat, "ichol_crash.mat");
    Eigen::IncompleteCholesky<float> fact;
    fact.compute(mat);
}

Steps to reproduce

Run the above code after decompressing the matrix market file ichol_crash.mat.gz in the working directory.

What is the current bug behavior?

Assertion:

/usr/local/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:412: Eigen::DenseCoeffsBase<Derived, 1>::Scalar& Eigen::DenseCoeffsBase<Derived, 1>::operator[](Eigen::Index) [with Derived = Eigen::Map<Eigen::Matrix<float, -1, 1>, 0, Eigen::Stride<0, 0> >; Eigen::DenseCoeffsBase<Derived, 1>::Scalar = float; Eigen::Index = long int]: Assertion `index >= 0 && index < size()' failed.

When not compiled in debug mode, this results in memory corruption (corrupted size vs. prev_size from glibc).

Stack:

#4  0x000055555555aec6 in Eigen::DenseCoeffsBase<Eigen::Map<Eigen::Matrix<float, -1, 1, 0, -1, 1>, 0, Eigen::Stride<0, 0> >, 1>::operator[] (this=0x7fffffffdbe0, index=1104202)
    at /usr/local/include/eigen3/Eigen/src/Core/DenseCoeffsBase.h:412
#5  0x0000555555558c29 in Eigen::IncompleteCholesky<float, 1, Eigen::AMDOrdering<int> >::factorize<Eigen::SparseMatrix<float, 0, int> > (this=0x7fffffffde00, mat=...)
    at /usr/local/include/eigen3/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h:254
#6  0x0000555555557d52 in Eigen::IncompleteCholesky<float, 1, Eigen::AMDOrdering<int> >::compute<Eigen::SparseMatrix<float, 0, int> > (this=0x7fffffffde00, mat=...)
    at /usr/local/include/eigen3/Eigen/src/IterativeLinearSolvers/IncompleteCholesky.h:146

What is the expected correct behavior?

No crash, hopefully ;)

Thanks!