function [L,D] = sdf(A)

% L * D * L' = A

[n,garbage] = size(A)

for i = 1:n
 for j = 1:n
  D(i,j) = 0;
 end
end


for k = 1:n
 L(k,k) = 1;
 tmp = 0;
 for v = 1:k-1
  tmp = tmp + D(v,v) * (L(k,v))^2;
 end
 D(k,k) = A(k,k) - tmp;
 for j = k+1:n
  L(k,j) = 0;
  tmp = 0;
  for v = 1:k-1
   tmp = tmp + L(j,v) * D(v,v) * L(k,v);
  end
  L(j,k) = (A(j,k) - tmp) / D(k,k);
 end
end
