function T = mvnfactor(sigma)     % originally [T,p] = 
%MVNFACTOR  Do Cholesky-like decomposition, allowing zero eigenvalues
%   SIGMA must be symmetric.  In general T is not square or triangular.

[U,D] = eig((sigma+sigma')/2);
D = diag(D);

tol = max(D) * length(D) * eps;
t = (D > tol);
D = D(t);
T = diag(sqrt(D)) * U(:,t)';