EDIT. 1. The Tommy L 's method is not better than the naive method.
Indeed, the complexity of the calculation of $(kI+A)^{-1}$ is $\approx n^3$ blocks (addition-multiplication).
About the complexity of $P(D+kI)^{-1}P^{-1}=QP^{-1}$ (when we know $P,D,P^{-1}$); the complexity of the calculation of $Q$ is $O(n^2)$ and the one of $QP^{-1}$ is $\approx n^3$ blocks as above. Note that one works with a fixed number of digits.
- When the real $k$ can take many (more than $n$) values, one idea is to do the following
The problem is equivalent to calculate the resolvent of $A$: $R(x)=(xI-A)^{-1}=\dfrac{Adjoint(xI-A)}{\det(xI-A)}=$
$\dfrac{P_0+\cdots+P_{n-1}x^{n-1}}{a_0+\cdots+a_{n-1}x^{n-1}+x^n}$.
STEP 1. Using the Leverrier iteration, we can calculate the matrices $P_i$ and the scalars $a_j$
$P_{n-1}:=I:a_{n-1}:=-Trace(A):$
for $k$ from $n-2$ by $-1$ to $0$ do
$P_k:=P_{k+1}A+a_{k+1}I:a_k:=-\dfrac{1}{n-k}Trace(P_kA):$
od:
During this step, we must calculate the exact values of the $P_i,a_j$ (assuming that $A$ is exactly known), which requires a very large number of digits. Then the complexity of the calculation is $O(n^4)$ (at least) but it's done only one time.
STEP 2. We put $x:=-k_1,-k_2,\cdots$. Unfortunately, the time of calculation -with a fixed number of significative digits- of $R(-k)$, is larger than the time of calculation of $(-kI-A)^{-1}$ !!!