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.)