Gaussian elimination

From Rosetta Code
Revision as of 08:37, 12 February 2012 by rosettacode>Pirho
Task
Gaussian elimination
You are encouraged to solve this task according to the task description, using any language you may know.

Problem: Solve Ax=b using Gaussian elimination then backwards substitution. A being an n by n matrix. Also, x and b are n by 1 vectors. To improve accuracy, please use partial pivoting and scaling.

MATLAB

<lang MATLAB> function [ x ] = GaussElim( A, b)

% Ensures A is n by n sz = size(A); if sz(1)~=sz(2)

   fprintf('A is not n by n\n');
   clear x;
   return;

end

n = sz(1);

% Ensures b is n x 1. if n~=sz(1)

   fprintf('b is not 1 by n.\n');
   return

end

x = zeros(n,1); aug = [A b]; tempmatrix = aug;

for i=2:sz(1)


   % Find maximum of row and divide by the maximum
   tempmatrix(1,:) = tempmatrix(1,:)/max(tempmatrix(1,:));
   
   % Finds the maximum in column
   temp = find(abs(tempmatrix) - max(abs(tempmatrix(:,1))));
   if length(temp)>2
       for j=1:length(temp)-1
           if j~=temp(j)
               maxi = j; %maxi = column number of maximum
               break;
           end
       end
   else % length(temp)==2
       maxi=1;
   end
   
   % Row swap if maxi is not 1
   if maxi~=1
       temp = tempmatrix(maxi,:);
       tempmatrix(maxi,:) = tempmatrix(1,:);
       tempmatrix(1,:) = temp;
   end    
   
   % Row reducing
   for j=2:length(tempmatrix)-1
       tempmatrix(j,:) = tempmatrix(j,:)-tempmatrix(j,1)/tempmatrix(1,1)*tempmatrix(1,:);
       if tempmatrix(j,j)==0 || isnan(tempmatrix(j,j)) || abs(tempmatrix(j,j))==Inf
           fprintf('Error: Matrix is singular.\n');
           clear x;
           return
       end
   end
   aug(i-1:end,i-1:end) = tempmatrix;
   
   % Decrease matrix size
   tempmatrix = tempmatrix(2:end,2:end);

end

% Backwards Substitution x(end) = aug(end,end)/aug(end,end-1); for i=n-1:-1:1

   x(i) = (aug(i,end)-dot(aug(i,1:end-1),x))/aug(i,i);

end

end </lang>