Eg
A = sparse.eye(5)
s = qdldl.Solver(A)
rhs = np.random.randn(5, 2)
s.solve(rhs)
The pure python implementation would be something like
[s.solve(rhs[:,i]) for i in range(rhs.shape[1])]
Most of the numpy/scipy solve commands support this so it might be nice.