The closest positive answers to your question that I could find is for sparse diagonal perturbations (see below).
With that said, I do not know of any algorithms for the general case, though there is a generalization of the technique you mentioned for scalar shifts from SPD matrices to all square matrices:
Given any square matrix $A$, there exists a Schur decomposition $A=U T U^H$, where $U$ is unitary and $T$ is upper triangular, and $A+\sigma I = U (T + \sigma I) U^H$ provides a Schur decomposition of $A + \sigma I$. Thus, your precomputation idea extends to all square matrices through the algorithm:
- Compute $[U,T]=\mathrm{schur}(A)$ in at most $\mathcal{O}(n^3)$ work.
- Solve each $(A+\sigma I) x = b$ via $x := U (T +\sigma I)^{-1} U^H b$ in $\mathcal{O}(n^2)$ work (the middle inversion is simply back substitution).
This line of reasoning reduces to the approach you mentioned when $A$ is SPD since the Schur decomposition reduces to an EVD for normal matrices, and the EVD coincides with the SVD for Hermitian positive-definite matrices.
Response to update:
Until I have a proof, which I do not, I refuse to claim that the answer is "no". However, I can give some insights as to why it's hard, as well as another subcase where the answer is yes.
The essential difficulty is that, even though the update is diagonal, it is still in general full rank, so the primary tool for updating an inverse, the Sherman-Morrison-Woodbury formula, does not appear to help. Even though the scalar shift case is also full rank, it is an extremely special case since it commutes with every matrix, as you mentioned.
With that said, if each $D$ was sparse, i.e., they each had $\mathcal{O}(1)$ nonzeros, then the Sherman-Morrison-Woodbury formula yields an $\mathcal{O}(n^2)$ solve with each pair $\{D,b\}$. For example, with a single nonzero at the $j$th diagonal entry, so that $D=\delta e_j e_j^H$:
$$
[A^{-1}+\delta e_j e_j^H]^{-1} = A^{-1} - \frac{\delta A^{-1} e_j e_j^H A^{-1}}{1+\delta (e_j^H A^{-1} e_j)},
$$
where $e_j$ is the $j$th standard basis vector.
Another update: I should mention that I tried the $A^{-1}$ preconditioner that @GeoffOxberry suggested on a few random SPD $1000 \times 1000$ matrices using PCG and, perhaps not surprisingly, it seems to greatly reduce the number of iterations when $||D||_2/||A||_2$ is small, but not when it is $\mathcal{O}(1)$ or greater.