tridiag_ldlt_sl.m

LU decomposition routines in matlab

function [b, iflag] = tridiag_ldlt_sl( d, e, b);
%
% TRIDIAG_LDLT_SL solves a tridiagonal linear system A x = b,
%    using the LDLT decomposition of A computed by TRIDIAG_LDLT
%
% Usage
%       [b, iflag] = tridiag_ldlt_sl( d, e, b);
%
% input:
%        d, e contain the L,D factors in compact form
%               of A as computed by tridiag_ldlt
%
%        b      right hand side vector
% output:
   
%
%        b  is overwritten with the solution if iflag= 0.
%
%        iflag: error flag
%               iflag = 0  solution could be computed and is
%                          stored in b
%               iflag = 1  dimension of b, and A are not compatible
%               iflag > 1  zero 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(b(:),1) )
   iflag = 1;
   return
end

% solve L y = b (overwrite b with the solution)
for i = 2:n
   b(i) = b(i) - e(i-1)*b(i-1);
end

% solve DL^T x = y (overwrite b with the solution)
b(n) = b(n) / d(n);
for i = n-1:-1:1
   b(i) = b(i)/d(i) - e(i)*b(i+1);
end