% 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