Let's say I have a sparse symmetric real indefinite matrix A of size N=150,000. This matrix is generated at each iteration of an algorithm with the same structure (positions of non-zero elements) but different values. I need to solve at each iteration several systems of the form Ax=b and to compute the matrix inertia (numbers of positive, negative and zero eigenvalues). From available literature, using an LDL.T factorization of the matrix A seems a promising way.
A is a KKT matrix, of the following form:
|--------+-----| | | | | H | G.T | | | | |--------+-----| | G | 0 | |--------+-----|
With H mostly diagonal, and G is rectangular, made of overlapping rectangular blocks along the diagonal (obtained by discretization of a state equations system).
I am using python. I have tried using the numpy package, the cvxopt package for dense and sparse implementations, but the factorization is way too slow for my goal.
Has anybody a python package / method to recommend for such a (medium) size problem? For example, a good implementation of a multifrontal direct solver/factorizer? Is there any way to exploit to the fullest the fact that the structure of the matrix is the same from one iteration to the next? Any thoughts outside of the box?