I would like to solve a series of linear systems Ax=b as quickly as quickly as possible. However, the systems are related. Specifically, each matrix A is given by:

cI + E

where E is a fixed sparse, symmetric positive definite real matrix (unchanged in all the linear systems), I is the identity matrix, and c is a varying complex number.

In other words, I am wondering how to quickly solve a series of complex linear systems which are all identical except for complex perturbations along the diagonal. I should say that the resulting matrices are not necessarily Hermitian, so currently I compute the LU decomposition. This works, but given the large number of rather closely related systems to be solved, I wonder if there is a better way to solve the problem, perhaps by using a more expensive (e.g. QR) decomposition up front.

(Edit for Jiahao: Yes, the bs are all the same.) (Edit for J. Mangaldan: The matrices are of order n=10^5 ~ 10^6, with about 10 times that many nonzeros.)

Update:

I'd like to thank everyone here for their suggestions. My implementation is ugly, but in the end interpolation was the key to a reasonable (10x) speedup. Since the c are quite close (imagine a small region of the complex plane, small in the sense that the spectrum of the matrix E is much larger) I could get away with computing solutions for a subset of the values of c and interpolating a solution for a given value of c using the precomputed values. It isn't elegant at all but it's something.