gauss.m

LU decomposition routines in matlab

function [A,b,iflag] = gauss( A, b)
%
% Usage:
%        [A,b,iflag] = gauss( A, b)
%
% Given an N by N matrix A and an N vector b,
% apply Gaussian elimination with column pivoting
% to transform the augmented matrix (A:b) into upper 
% triangular form
%
% input:
%        A:     the N by N matrix A
%
%        b:     the N vector b
%
% output:
%        A:     the transformed N by N matrix A
%
%        b:     the transformed N vector b (matrix b)
%
%        iflag: error flag
%               iflag = 0  Row reductions could be performed,
%                          A is upper triangular
%               iflag = 1  dimension of A or of b is not correct
%
%
% Matthias Heinkenschloss
% Sept. 22, 1995

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 Gaussian elimination

for k = 1:n-1

%   find row index containing the element in column
%   k that has largest absolut value

    [amax,i] = max(abs(A(k:n,k)));
    i = i + k - 1;
 
%   perform row interchange if necessary
    if( i ~= k ) 
        atmp = A(k,k:n);
        A(k,k:n) = A(i,k:n);
        A(i,k:n) = atmp;
        btmp   = b(k,:);
        b(k,:) = b(i,:);
        b(i,:) = btmp;
    end 

    if( A(k,k) ~= 0)
%       perform row elimination
        for i = k+1:n    
%           compute multiplier
            l = A(i,k)/A(k,k);

            A(i,k:n) = A(i,k:n) - l * A(k,k:n);
 
            b(i,:) = b(i,:) - l * b(k,:);
            A(i,k) = 0;
        end
    end
end