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.
Best Answer
You want the resolvent of $E$ (at $z:=-c$). Recall it's an analytic function of $z$ defined on the resolvent set, $\mathbb{C}\setminus\operatorname{spec}(E)$. According to the complexity of the matrix $E$, and with the number and the location of the $c$ you need to consider, it may be worth computing a power series expansion at various centers so as to cover the set $\{c\}$ of the data. For $|z|$ larger than the spectral radius you have of course the Laurent expansion $(z-E)^{-1}=1/z+ E/z^2+ E^2/z^3+\dots$