function [A,D] = triang_fact(Omega);
%------------------------------------------------------------
%Purpose: calculates the unique triangular factorization 
%         of a real symmetric positive definite matrix Omega
%------------------------------------------------------------     
%      Omega = A*D*A'    A is lower triangular with ones down
%                        the principal diagonal
%                        D is a diagonal matrix with positive 
%                        entries down the main diagonal
%------------------------------------------------------------
% Written by: Gregory Givens
%             3 March 2003
%------------------------------------------------------------

%check to see if Omega is square symmetric positive definite matrix
[m,n] = size(Omega); F = Omega';
if m ~= n;
    error('matrix is not square');
end;
for i = 1:m;
    for j = 1:n;
        if Omega(i,j) ~= F(i,j);
            error('matrix is not symmetric');
        end;
    end;
end;
Lambda = eig(Omega);
for k = 1:m;
    if Lambda(k) <= 0;
        error('matrix is not positive definite');
    end;
end;

P = chol(Omega)'; D_s = diag(diag(P));
A = P*inv(D_s); D = D_s*D_s;