function [L,D] = ldldecom(Q)
%LDLDECOM: Find LtDL-decompostion of Q-matrix
%
% This routine finds the LtDL decomposition of a given variance/
% covariance matrix.
%
% Input arguments:
%    Q: Symmetric n by n matrix to be factored
%
% Output arguments:
%    L: Out - n by n factor matrix (strict lower triangular)
%    D: Out - Diagonal n-vector

% ----------------------------------------------------------------------
% File.....: ldldecom
% Date.....: 19-MAY-1999
% Author...: Peter Joosten
%            Mathematical Geodesy and Positioning
%            Delft University of Technology
% ----------------------------------------------------------------------

n = size (Q,1);

for i = n:-1:1;

   D(i) = Q(i,i);
   L(i,1:i) = Q(i,1:i)/sqrt(Q(i,i));
   
   for j = 1:i-1
      Q(j,1:j) = Q(j,1:j)-L(i,1:j)*L(i,j);
   end
   
   L(i,1:i) = L(i,1:i)/L(i,i);

end;

if (sum(D < 1E-10));

  error ('Matrix on input is not positive definite!');

end;

% ----------------------------------------------------------------------
% End of routine: ldldecom
% ----------------------------------------------------------------------