function [A, r] = gauss(A)
% gauss Scaled maximal column pivoting with Gaussian elimination.
% right hand side vector is not included
% Syntax: [A,r] = gauss(A)
% A is the n by n matrix (changed by gauss)
% r is the pivot vector returned for use with gaussSolve

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

% initialize the pivot vector

r = zeros(1,n);
for i = 1 : n
   r(i) = i;
end

% First compute scale factors scale[1] to scale[n]
% for each row

scale = zeros(1,n);
for i = 1 : n
   smax = 0;
   for j = 1 : n
      w = abs(A(i,j));
      if w > smax
         smax = w;
      end
   end
   if smax == 0
      error('matrix is singular');
   end
   scale(i) = smax;
end

% forward elimination

% loop over pivot positions (diagonal elements)
for k = 1 : n-1
   
   % to index matrix element A(i,j) we now need to
   % use A(r(i),j) since rows are not exchanged
  
   p = k;
   amax = 0;
   for i = k : n
      w = abs(A(r(i),k)/scale(r(i)));
      if w > amax;
         amax = w;
         p = i;
      end
   end
   
   % if we didn't find a non-zero pivot then
   % the matrix is singular
   
   if amax == 0
      error('Matrix is singular');
   end
   
   % if we found a non-zero pivot in row p then exchange
   % pivot vector entries for row k and row p
   
   if p > k
      temp = r(k);
      r(k) = r(p);
      r(p) = temp;
   end

   % Note: we don't need to exchange elements in the
   % scale vector because they will be accessed using
   % scale(r(i)) instead of scale(i)
   
   % loop down pivot column k beginning below pivot
   for i = k+1 : n
      
      % calculate multiplier for row r(i)
      % and subtract multiple of row r(i) from row r(k)
      % Note: we store the mutipliers in the lower
      % triangular part of the matrix
      
      mult = A(r(i),k) / A(r(k),k);
      for j = k+1 : n
         A(r(i),j) = A(r(i),j) - mult*A(r(k),j);
      end
      A(r(i),k) = mult;
   end
end

end % gauss