Let $\Sigma$ be positive definite and the $D_i$ positive diagonal. Let the $X_i$ be unknown square matrices. Consider the system of equations:

$$(I+\Sigma D_i)X_i=\Sigma\hspace{5mm}\text{for}\hspace{5mm}i=1,...n.$$

Is it possible to numerically solve all of these equations without having to completely recompute the LU factorization of $(I+\Sigma D_i)$ for each $i$?

I've been messing around with the Cholesky and Eigen-decompositions of $\Sigma$, but no luck so far, I think it's probably not possible, but I thought I'd ask.

If it helps, note that the solutions $X_i$ will all be positive definite, since

$$X_i^{-1} = \Sigma^{-1} + D_i.$$

Edit:

It looks like this answer here by the all knowing Brian Borchers means that this is infact impossible.