Implementation of Partial Computed Gradients (PCG) Method for Distributed, Sparse Arrays

I'm running a transient FEM of a very large matrix (>1E6x1E6) on a cluster, and looking for a way to solve it the most efficient way possible. The matrix must be sparse, because a full matrix would exceed memory. Is there a way to use distributed arrays on this sparse data to solve the matrix, using the pcg method? Using just sparse and pcg, each time step of the simulation takes 10-30 seconds, and I'm looking at >200,000 time steps. If I can harness multiple cores to speed this up by even a factor of 2, then I would be very happy.