Sorry about the long title. I need to calculate the trace of $I(I+P)^{-1}$$M(M+D)^{-1}$, where $I$$M$ is a dense symmetric matrix, and $P$$D$ is a diagonal matrix. The main issue is the dimension could be large (usually in the hundreds to thousands range) so computing the inverse could be quite expensive especially because the trace needs to be computed iteratively. Is there any efficient algorithm to do this?