Is there some sort of memory-efficient approximate solution? What if all of the off-diagonal elements of the covariance matrix are equal to $\alpha$ and all of the diagonal elements are equal to $\sigma^2$? An implementation in R would be terrific.
For additional context, the MVN has a dimension around 70000 and has a dense covariance matrix. I've tried using the standard multivariate normal packages in R but even on a machine with a lot of memory can't realistically sample from MVNs with more than a few thousand elements. I'm running simulations to test for coverage and power of parameter estimates in settings where there's high covariance between observations, and it would be helpful if I could exactly specify the covariance matrix I want to use.