function x = geppbs(A,b,ppon)
% GEPPBS  Gaussian Elimination with partial pivoting and backward substitution.
% Reasonable implementation for serious work, but of course it duplicates "A\b"
% entirely.  Note that actually moving rows is unnecessary.
%
% Example 1; basic:
%   >> A = [3 2 -1; -1 1 1; 2 0 2]          % enter system of equations
%   >> b = [4 -1 5]'
%   >> x = geppbs(A,b)
%   >> X = A \ b                            % compare to built-in
% Example 2; show effect of partial pivoting:
%   >> A = randn(500,500);  x = randn(500,1);     % random system ...
%   >> b = A * x;                                 %    known solution
%   >> X = geppbs(A,b,true);  norm(x-X)/norm(x)   % relative error with p.p.
%   >> X = geppbs(A,b,false);  norm(x-X)/norm(x)  % relative error without p.p.
%
% See also:  draft_ge.m, gebueler.m, geexample.m

if nargin < 3, ppon = true; end  %  if third argument not given then do p.p.

n = size(A,1);              % number of rows = number of equations
if size(A,2) ~= n           % check on reasonable shape of augmented matrix
  error('A must be square (n x n)')
end
if (size(b,1) ~= n) | (size(b,2) ~= 1)
  error('b must be  n x 1  if A is  n x n')
end

% GE = Gaussian elimination with partial pivoting phase
for i = 1:n-1                      % go through cols in which to make zeros
  p = i;
  if (ppon)
    for k = i+1:n                  % partial pivoting finds  p  so that
      if abs(A(k,i)) > abs(A(p,i)) %     |A(p,i)| >= |A(k,i)|  for  k = i:n
        p = k;
      end
    end
  else
    while A(p,i) == 0.0
      p = p+1;
    end
  end
  if A(p,i) == 0.0, error('no soln'), end % fail if everything in col was zero
  if p > i                         % do swap if needed
    tmp = A(i,:);
    A(i,:) = A(p,:);
    A(p,:) = tmp;
    tmpb = b(i);
    b(i) = b(p);
    b(p) = tmpb;
  end
  for j = i+1:n                    % go through rows below diagonal
    m = A(j,i) / A(i,i);           % multiplier
    A(j,:) = A(j,:) - m * A(i,:);  % operation:   E_j <- E_j - # E_i
    b(j) = b(j) - m * b(i);        % same action on RHS
  end
end 
if A(n,n) == 0, error('no soln'), end

% BS = backward substitution phase
x = zeros(size(b));                % allocate the right space
x(n) = b(n) / A(n,n);
for i = n-1:-1:1
  x(i) = (b(i) - A(i,i+1:n) * x(i+1:n)) / A(i,i);
end

