ltriangsl.m

LU decomposition routines in matlab

function [b,iflag] = ltriangsl( A, b);
%
% LTRIANGSL  solves an N by N linear system with
% lower triangular matrix.
%
% Usage
%       [b,iflag] = ltriangsl( A, b)
%
% input:
%        A:     the lower triangular N by N matrix A
%
%        b:     the right hand side b. b can have more than one column. 
%               In this case a system with multiple right hand sides is 
%               solved.
%
%
% output:
%        b:     the solution of the system A x = b, if iflag = 0
%
%        iflag: error flag
%               iflag = 0  solution could be computed and is
%                          stored in b
%               iflag = 1  dimension of A or of b is not correct
%               iflag > 1  zero diagonal element detected in row iflag+1
%
%


iflag = 0;

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

% Start forward substitution (column form)

if( A(1,1) == 0 )
    iflag = 2;
    return
end
b(1,:) = b(1,:) / A(1,1);

for i = 1:n-1
   b(i+1:n,:) = b(i+1:n,:) - A(i+1:n,i) * b(i,:);
   if( A(i+1,i+1) == 0 )
       iflag = i+2;
       return
   end 
   b(i+1,:) = b(i+1,:) / A(i+1,i+1);
end