  function Show_LDL_LDM()
% Illustrates the A= LDL' and A= LDM' factororizations.
% To ensure stability we assume that A is sym pos def in LDL'
% and diagonally dominant in LDM'. 

n = 5;
% Generate a small random symmetric positive definite example
A = randn(n,n); A = A'*A
[L,d] = LDL(A)
errLDL = norm(A-L*diag(d)*L');
% Generate a small random diagonally dominant example
A = randn(n,n); 
A = A + 2*diag(sum(abs(A)))
[L,d,M] = LDM(A)
errLDM = norm(A-L*diag(d)*M');

fprintf('errLDL = %10.2e\nerrLDM = %10.2e\n',errLDL,errLDM)

  function [L,d] = LDL(A)
% A is an nxn and symmetric positive definite matrix
% L is unit lower triangular and d is a column n-vector such that
% A = LDL' where D = diag(d).
[n,n] = size(A);
d = zeros(n,1);
L = eye(n,n);
for j=1:n
    % Compute d(k) and the k-th column of L...
    s          = d(1:j-1).*L(j,1:j-1)';
    d(j)       = A(j,j) - L(j,1:j-1)*s;
    v(j+1:n)   = A(j+1:n,j) - L(j+1:n,1:j-1)*s;
    L(j+1:n,j) = v(j+1:n)/d(j);
end

  function [L,d,M] = LDM(A)
% A is nxn and has an LDM factorization
% L and M are nxn and unit lower triangular and d is a column n vector 
% such that A = LDM' where D = diag(d).
[n,n] = size(A);
d = zeros(n,1);
L = eye(n,n);
M = eye(n,n);
for j=1:n
    % Compute d(k), the kth column of L and the kth column of M...
    s    = d(1:j-1).*M(j,1:j-1)';
    d(j) = A(j,j) - L(j,1:j-1)*s;
    v(j+1:n) = A(j+1:n,j) - L(j+1:n,1:j-1)*s;
    L(j+1:n,j) = v(j+1:n)/d(j);
    t    = d(1:j-1).*L(j,1:j-1)';
    w(j+1:n) = A(j,j+1:n) - t'*M(j+1:n,1:j-1)';
    M(j+1:n,j) = w(j+1:n)/d(j);
end
