% gauss_sp.m: Gaussian elimination with scaled partial pivoting.

clear;
format short;

% Step 0: Assign the matrix A and the vector b.

n = 4;
A = [ 6,  -2, 2,   4;
12,  -8, 6,  10;
3, -13, 9,   3;
-6,   4, 1, -18 ];
b = [ 0.23; 0.32; 0.33; 0.31 ];

% Step 1: Gaussian Elimination with scaled partial pivoting.
%         Overwrite the matrix.

p = (1:n)';	        % initialize the pivoting vector
s = max(abs(A'));     % compute the scale of each row
for k = 1:(n-1)
r = abs(A(p(k),k)/s(p(k)));
kp = k;
for i = (k+1):n
t = abs(A(p(i),k)/s(p(i)));
if t > r,  r = t;  kp = i;  end
end
l = p(kp);  p(kp) = p(k);  p(k) = l;   % interchange p(kp) and p(k)
for i = (k+1):n
A(p(i),k) = A(p(i),k)/A(p(k),k);
for j = (k+1):n
A(p(i),j) = A(p(i),j)-A(p(i),k)*A(p(k),j);
end
end
end
p                      % output the pivoting vector
A                      % output the overwritten matrix

% Step 2: Forward subsitiution to solve  L*y = b, where
%         L(i,j) = 0 for j>i; L(i,i) = 1;  L(i,j) = A(p(i),j) for i>j.

y = zeros(n,1);        % initialize y to be a column vector
y(1) = b(p(1));
for i = 2:n
y(i) = b(p(i));
for j = 1:(i-1)
y(i) = y(i)-A(p(i),j)*y(j);
end
end
y                      % output y

% Step 3: Back subsitiution to solve  U*x = y
%         U(i,j) = A(p(i),j) for j>=i; U(i,j) = 0 for i>j.

x = zeros(n,1);        % initialize x to be a column vector
x(n) = y(n)/A(p(n),n);
for i = (n-1):-1:1
x(i) = y(i);
for j = (i+1):n
x(i) = x(i)-A(p(i),j)*x(j);
end
x(i) = x(i)/A(p(i),i);
end
x                      % output x