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