tridiag_lu_sl.m

LU decomposition routines in matlab

function [b, iflag] = tridiag_lu_sl( c, d, e, b);
%
% TRIDIAG_LU_SL solves a tridiagonal linear system A x = b,
%    using the LU decomposition of A computed by TRIDIAG_LU.
%
% Usage
%       [b, iflag] = tridiag_lu_sl( c, d, e, b)
%
% input:
%        d, e contain the LU factors in compact form
%               of A as computed by tridiag_lu
%
%        b      right hand side vector
% output:
   
%
%        b  overwritten with the solution of the system 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) - c(i-1)*b(i-1);
end

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