function [U,D]=udut(R)
% function [U,D]=udut(R)
% Computes the Upper-Lower UDL 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);
U=zeros(M,M);
for m=1:M U(m,m)=1; end
D=zeros(M,1);

D(M)=R(M,M);

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