ldlt.m

LU decomposition routines in matlab

function [A,iflag] = ldlt( A )
%
% Usage:
%        [A,iflag] = ldlt( A )
%
% Given a symmetric N by N matrix A, LDLT attempts
% to compute the LDL^T decomposition of A.
% 
%
% input:
%        A:     The N by N matrix A.
%               Only the lower triangular part of A
%               is used.
%
% output:
%        A:     The strict lower triangular part of A
%               contains the strict lower triangular part of L.
%               The diagonal of A contains the diagonal D.
%               The strict upper part of A is unchanged.
%
%        iflag: error flag
%               iflag = 0  LDL^T decomposition could be computed
%                          and is stored in the lower triangular
%                          part of A.
%               iflag = 1  A is not squared
%               iflag > 1  Negative diagonal entry observed in step iflag-1.
%                          Matrix A is not positive definite.
%
%
% Matthias Heinjenschloss
% Feb 1, 2001

iflag = 0;

% get size of A and checj dimensions
[m,n]    = size(A);
if ( m ~= n )
   iflag = 1;
   return
end


% Start LDL^T decomposition

for j = 2:n

    A(j,j) = A(j,j) - A(j,1:j-1)*A(j,1:j-1)';
    if( A(j,j) <= 0 )
        iflag = j+1;
        return
    end
    
    for i = j+1:n
        A(i,j) = ( A(i,j) - A(i,1:j-1)*A(j,1:j-1).*' ) / A(j,j);
    end
end
    
%end of ldlt