RosettaCodeData/Task/Gaussian-elimination/MATLAB/gaussian-elimination.m

71 lines
1.5 KiB
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