function [L,D]=ldlt(R)
% function [L,D]=ldlt(R)
% Computes the LDU decomposition of a 
% symmetric positive definite matrix.
%
% Programmed by: Dimitris Manolakis, 1996
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------



M=length(R);
L=zeros(M,M);
for m=1:M L(m,m)=1; end
D=zeros(M,1);

D(1)=R(1,1);
L(2,1)=R(2,1)/D(1);
D(2)=R(2,2)-D(1)*L(2,1)*L(2,1);

for j=3:M
  for i=1:j-1
    L(j,i)=R(j,i);
    for k=1:i-1 L(j,i)=L(j,i)-D(k)*L(j,k)*L(i,k); end
    L(j,i)=L(j,i)/D(i);
  end
  D(j)=R(j,j);
  for k=1:j-1 D(j)=D(j)-D(k)*L(j,k)*L(j,k); end
end