function [x] = zgauss(A, b)
% zgauss gaussian elimination.
% This version tries to avoid zero pivots
% Syntax: x = zgauss(A,b)
% A is the n by n matrix (not changed by zgauss)
% b is the right hand side vector (not changed by zgauss)
% x is the solution (column vector)

[m,n] = size(A);
if m ~= n
   error('Matrix is not square');
end

% forward elimination

% loop over pivot positions (diagonal elements)
for k = 1 : n-1
   
   % if pivot element is zero try to find a non-zero
   % entry below it in column k
   
   p = k;
   while p <= n && A(p,k) == 0
      p = p + 1;
   end
   
   % if we didn't find a non-zero pivot then
   % the matrix is singular
   
   if p > n
      error('Matrix is singular');
   end
   
   % if we found a non-zero pivot in row p then exchange
   % pivot row k with row p
   
   if p > k
      for j = 1 : n
         temp = A(k,j);
         A(k,j) = A(p,j);
         A(p,j) = temp;
      end
      temp = b(k);
      b(k) = b(p);
      b(p) = temp;
   end
   
   
   % loop down pivot column k beginning below pivot
   for i = k+1 : n
      
      % calculate multiplier for row i
      % and subtract multiple of row k from row i
      mult = A(i,k) / A(k,k);
      for j = 1 : n
         A(i,j) = A(i,j) - mult*A(k,j);
      end
      
      % do the same for the right hand side vector
      b(i) = b(i) - mult*b(k);
   end
end

% Back substitution beginning with last unknown x(n)
% and proceding backward until x(1) is found

x = zeros(n,1);
x(n) = b(n) / A(n,n);
for i = n-1 : -1 : 1;
   % compute unknown x(i) in terms of x(i+1) to x(n)
   s = b(i);
   for j = i+1 : n
      s = s - A(i,j)*x(j);
   end
   x(i) = s / A(i,i);
end

end % zgauss