function  [X_pd, iflag] = make_pd(X)
    small = 1.0e-8;
    % Check if X is positive definite
    [~,p] = chol(X);
    if p == 0
        X_pd = X;
        iflag = 0;
    else
        iflag = 1;
        scl = mean(diag(X));
        Y = X/scl;
        % Matrix is not PD ... make it PD by truncating small eigenvalues
        [V,D] = eig(Y);
        d = diag(D);
        ii = d<small;
        d(ii==1) = small;
        D1 = diag(d);
        V = real(V);
        Y_pd = V*D1*V';
        X_pd = Y_pd*scl;
    end
     
end