tridiag_ldlt.m

LU decomposition routines in matlab

function [d, e, iflag] = tridiag_ldlt( d, e);
%
% TRIDIAG_LDLT computes the LDL^T-decomposition of a 
%    symmetric positive definite tridiagonal matrix.
%
% Usage
%       [c, d, e, iflag] = tridiag_lu( c, d, e)
%
% input:
%        c:     real(n-1)
%               the subdiagonal entries of A 
%
%        d:     real(n)
%               the diagonal entries of A 
%
%        e:     real(n-1)
%               the superdiagonal entries of A 
%
%
%        The matrix is stored as
%
%               d(1) e(1)
%               e(1) d(2) e(2)
%                    e(2) d(3) e(3)
%                       .    .    .
% output:
   
%
%        d, e are overwritten with the L D factors in compact form
%               of A
%
%
%        iflag: error flag
%               iflag = 0  LU decomposition could be computed and is
%                          stored in  d, e
%               iflag = 1  dimension of d, e are not correct
%               iflag > 1  zero or negative diagonal element detected 
%                          in step iflag+1
%
%
%
% Matthias Heinkenschloss
% Feb 6, 2001


iflag = 0;

% get size of A and check dimensions
n    = size(d(:),1);
if ( n > size(e(:),1)+1 )
   iflag = 1;
   return
end

for i = 1:n-1
   if( d(i) <= 0 )
       iflag = i+1;
       return
   end 
   t = e(i);
   e(i)   = e(i) / d(i);
   d(i+1) = d(i+1) - t*e(i);
end