function r=sample_multinorm(mu,sigma)
% Copyright (c) 2013-2018, Francesco Corsello, All Rights Reserved.
[n,d] = size(mu);
sz = size(sigma);
if n>1 || not(ismatrix(sigma)) || ~isequal(sz, [d d])
    error('Mu must be a row vector and Sigma must be a matrix of the correct size.')
end     

        [T,p] = chol(sigma);
        if p>0            
            [U,D] = eig((sigma+sigma'/2));
            [~,maxind] = max(abs(U),[],1);
            negloc = (U(maxind + (0:n:(n-1)*n)) < 0);
            U(:,negloc) = -U(:,negloc);
            D = diag(D);
            tol = eps(max(D)) * length(D);
            t = (abs(D) > tol);
            D = D(t);            
            T = diag(sqrt(D)) * U(:,t)';
        end
        r = randn(n,size(T,1)) * T + mu;