Solving large-scale symmetric linear system with matrix inertia (python)

My requirements are to solve the following linear system and determine the LHS matrix inertia using python:

(Matrix Inertia = (p,n,z) - The number of positive, negative and zero eigenvalues of the RHS matrix.)

Large-Scale Linear System

The solution to this system derives the search direction for a primal-dual interior point optimisation.

Each term in this figure has the following properties:

n >> 1000
m < 10
W = [n x n] (Sparse)
Sigma = [n x n] (Sparse)
delta_w_I = [n x n] (Sparse)
A = [n x m]
delta_c = [m x m]
dx = [n x 1]
d_lambda = [m x 1]
delta_phi = [n x 1]
A_lambda = [n x 1]
c = [m x 1]

W represents the hessian and must be approximated using the L-BFGS which only computes the inverse W^-1.

What steps do I need to take to solve this sparse system efficiently in python whilst also calculating the matrix inertia?

I have read about LDL^T factorisation but it doesn't work as I do not wish to invert my W term as this method requires both W and W^-1 to solve and determine the matrix inertia.

(I have tried to install the Harwell subroutine libraries on my machine with the HSL.py wrapper, but I find it impossible to do so.)