Symmetric Doolittle Factorization

sdf.m

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

L, D, L * D * L'

L =

1.0000         0         0
0.5000    1.0000         0
     0    2.0000    1.0000

D =

 4     0     0
 0     2     0
 0     0     0

ans =

 4     2     0
 2     3     4
 0     4     8