Suppose $$\bf{\nu}|\mu,\rho,\delta^2 \sim N_p(\mu \bf{1},\delta^2(D_w-\rho W)^{-1}),$$ where $W$ is a binary symmetric matrix and $D_w$ is diagonal with $(D_w)_{ii}=\sum_j w_{ij}$, $\mu$ is a scalar.
How to prove that
$$cov(\nu_i,\nu_j|\rho,\delta^2,{\{\nu_k\}}_{k\neq i,j})=\frac{\delta^2 \rho w_{ij}}{w_{i0} w_{j0}+\rho^2 w_{ij}^2},$$ where $w_{i0}$ and $w_{j0}$ are the $i^{th}$ and $j^{th}$ diagonal elements of $D_w$ respectively ?
[As a reference see Gelfand and Vounatsou 2003 article.]