Good evening, i've written a program that uses Eigen3 to solve a sparse linear system where the input matrix is an SPD matrix in .mtx format and the output x should be a vector of ones. I've to test 9 different matrices, the program works fine with the first 7 matrices, while the 8th matrix causes a "write access violation" exception. The first 7 matrices are all < 100MB of dimension, while this one is around 300MB.
This is the code:
typedef Eigen::SimplicialLDLT<SM> CS;
typedef Eigen::VectorXd V;
typedef Eigen::SparseMatrix<double> SM;
int
main (int argc, char *argv[])
{
SM mat;
Eigen::loadMarket(mat, std::string(argv[1]));
SM A = mat.selfadjointView<Eigen::Lower>();
CS solver;
V b(A.rows(), 1), x(A.rows(), 1), xe(A.rows(), 1);
xe.setOnes(A.cols(), 1);
b = A * xe;
solver.compute(A);
x = solver.solve(b);
}
The crash occours on solver.compute(A). Debugging the code i've found out that the error is inside SimplicialCholesky_impl.h
Li[p] = k; /* store L(k,i) in column form of L */
Li is defined as follows:
StorageIndex* Li = m_matrix.innerIndexPtr();
but the value of Li is 0x0000000000000000 and i suppose this is somehow wrong, but i can't really understand what to do to solve the problem, especially because this happens only with this specific matrix.
The interested matrix is this one https://www.cise.ufl.edu/research/sparse/matrices/Janna/StocF-1465.html .
I'm working with msvc and visual studio 2017, the build is release x64.
By default SparseMatrix uses int to store indices and thus so does SimplicialLDLT<SM> for its L factor. For your problem, you clearly need long int, so all you have to do is:
typedef Eigen::SparseMatrix<double,ColMajor,long> SM;
but this will take time because non-supernodal Cholesky factorization are good for 2D problems only and this matrix comes from a 3D finite element discretization.