function [L,d,p,e] = mchol(H)
% Modified Cholesky Factorization (see GMW pg 111).
% P'H P + E = L D L'
% L is lower triangular, D = diag(d), E = diag(e)
% p is permutation vector p_j = i iff P_ij = 1

delta = 1e-8;
[n,m] = size(H);
if (n ~= m)
    error('Cholesky needs square matrices');
end
nu = max(1.0,sqrt(n*n - 1.0));
gamma = max(abs(diag(H)));
eta = max(abs(triu(H,1)));
beta_2 = max([gamma,eta/nu,eps]);
d = zeros(n,1);
e = zeros(n,1);
p = 1:n;
C = diag(diag(H));
L = speye(n);
for j=1:n,
    [cqq,q] = max(abs(diag(C)));
    if (q>j) 
      %  swap rows and cols q and j of H and C
      p([j q]) = p([q j]);
      C([j, q],:) = C([q,j],:);
      C(:,[j, q]) = C(:,[q,j]);
    end

    k = 1:j-1;
    i=j+1:n;
    if isempty(k)
      C(i,j) = H(p(i),p(j));
    else
      L(j,k) = C(j,k)./d(k)';
      C(i,j) = H(p(i),p(j)) - C(i,k)*L(j,k)';
    end
    theta_j = max([0;abs(C(i,j))]);
    d(j) = max([delta,abs(C(j,j)),theta_j*theta_j/beta_2]);
    e(j) = d(j) - C(j,j);
    C(i,i) = C(i,i) - C(i,j)*C(i,j)'/d(j);
    C(j,j) = 0.0;
end
return;