Gaussian Elimination Script

Question

I'm trying to use the following script:

> function [x,An,Ln,counter]=gausselim(A,b,details)
% This function performs Gaussian elimination
% on the square matrix A with right hand side b.
%
% INPUT: square (n,n) matrix A
%        (n,dimb) vector b
%         details: 1 if intermediate steps should be shown
%                  0 if not.
% OUTPUT solution vector x of equation Ax=b;
%        upper triangular matrix An after n-1 steps
%        lower triangular matrix Ln containing mjl
%        number of row exchanges

% Notice that the function written here is not very
% efficient with respect to saving memory: Typically
% we would not have to define extra matrices Ln and An
% for this procedure.

n=size(A,1); 
if size(A,2)~=n 
    error('input matrix not square');
end
dimb=size(b,2);
if size(b,1)~=n 
    error('wrong dimension of right hand side');
end
% Form the augmented matrix Atilde=[A,b]
% and call it A
A=[A,b];
% print matrix if details are required
fprintf('\n');
if details==1
       fprintf('This is the matrix A%d: \n', 1); 
       A
end
% In Ln the multiplier mjl are stored (and 1 on diagonal)
Ln=eye(n);
% Initialize solution vectors x 
x=zeros(n,dimb);
% Initialize counter which counts row exchanges:
counter=0;
% Step G1:
% Perform elimination steps EL1-EL6:
for l=1:(n-1) 
% Step EL1: look for nonzero pivot element in column l
    p=-1;
    j=l;
    while p==-1 
       if A(j,l)~=0 
           p=j;
       else 
           j=j+1;
       end    
       if j>=(n+1)
       fprintf('error: pivot is zero in step %d\n',l);
       fprintf('No unique solution can be found\n');
       An=A(:,1:n);
         return  
       end
    end
% Step EL2: exchange rows if necessary
   if p~=l
   counter=counter+1;
   pivotrow=A(p,:);
   A(p,:)=A(l,:);
   A(l,:)=pivotrow;
   end
% Step EL3: eliminate elements in column l below Pivot:
   for j=(l+1):n
   mjl=A(j,l)/A(l,l);
   Ln(j,l)=mjl;
   A(j,:)=A(j,:)-mjl.*A(l,:);
   end
% Step EL6: nothing is to be done here
  ;

end     
% Step G2: Check whether n-th pivot is nonzero
if A(n,n)==0
       fprintf('error: pivot is zero in step %d\n',n);
       fprintf('No unique solution can be found\n');
       An=A(:,1:n);
        return
end      
% This completes the elimination step.
% Write upper triangular matrix An
An=A(:,1:n);

% Backward substitution:
% perform on each vector b separately:

% Step B3:
disp('Gaussian elimination completed successfully');
fprintf('it needed %d row exchanges\n\n', counter);

But no matter what matrix I try to use it with I get something like the following:

Gaussian elimination completed successfully
it needed 0 row exchanges


x =

     0
     0
     0
     0
     0

Test: the product A*x for this x is:

ans =

     0
     0
     0
     0
     0

Can someone please tell me what I'm doing wrong?

Thanks.

Answer 1

Well, you never update your vector x. You initialize your x to x=zeros(n,dimb); in line 39 and you never update that. Obviously you will get x = [0, ..., 0] at the end of the algorithm.