function [x, m, L] = fastmvg_rue(Phi, PtP, alpha, D, XtX, gprior)
%FASTMVG_RUE sampler for multivariate Gaussians. 
%   [x, m, L] = fastmvg_rue(...) generates multivariate Gaussian random
%   variates of the form: N(mu, S), where
%       mu = S Phi' y
%       S  = inv(Phi'Phi + inv(D))
%
%   Here, PtP = Phi'*Phi (X'X may be precomputed). Sampler used for n
%   larger than p.
%   
%   The input arguments are:
%       Phi    - [n x p] matrix 
%       PtP    - [p x p] matrix
%       alpha  - [p x 1] vector
%       D      - [p x 1] vector
%       XtX    - [p x p] pre-computed X'*X (if available)
%       gprior - [1 x 1] true for gprior, otherwise false
%
%   Return values:
%       x      - [p x 1] vector of MVG random variates
%       m      - [p x 1] vector, posterior mean of the MVG 
%       L      - [p x p] Cholesky factor 
%
%   Reference:
%     Rue, H. (2001). Fast sampling of gaussian markov random fields. Journal of the Royal
%     Statistical Society: Series B (Statistical Methodology) 63, 325338.
%
%   (c) Copyright Enes Makalic and Daniel F. Schmidt, 2016

if(isempty(PtP))
    PtP = Phi' * Phi;
end

p = size(Phi,2);

if(~gprior)
    L = chol(PtP + diag(1./D), 'lower');
else
    L = chol(PtP + XtX./D(1), 'lower');
end
v = L \ (Phi'*alpha);
m = L' \ v;
w = L' \ randn(p,1);

x = m + w;

end