function y = randvec(m, mu, sigma)
%
% RANDVEC: y = randvec(m, mu, sigma);
%
% This routine returns m random vectors (one in each column of y).
% The vectors are distributed according to a multivariate normal 
% distribution with mean mu and covariance sigma. 
%
% Note: mu must be a column vector.
%

[p, r] = size(sigma);
x = randn(p, m);
%%q = (chol(sigma))';
[U,S,V] = svd(sigma);
q = U * sqrt(S);        % Changed from chol to svd to handle non-neg dfeinite
y = q * x + mu * ones(1,m);