function [A,b] = gaussianelimination(A,b)
% This function converts a linear system of equations with coefficient matrix
% $\VA$ and right hand side vector $\Vb$ into canonical for
% (``Zeilenstufenform'') by performing row transformations. 
% \texttt{A} must pass an $n\times m$-matrix and \texttt{b} a column vector of
% length $n$.
[n,m] = size(A); % Obtain number of rows $n$ and number of columns $m$ of $\VA$
% Remember that testing for equality is not a  good idea, if floating point
% operations are involved, because roundoff will render this test meaningless. 
% Therefore it is recommended that one tests relative equality up to a relative
% perturbation below a prescribed tolerance. This tolerance is specified in the
% next line. 
tol = eps*n*m;

% First trick: add $\Vb$ as rightmost column to $\VA$ in order to facilitate
% implementation of simultaneous row transformations of $\VA$ and $\Vb$
A = [A,b], 

% \texttt{j} is the current column index, \texttt{l} the current row index.
j = 1; l = 1;
while (j<=m)
  % Search from top of column \texttt{j} to find first non-zero entry.
  % In the absence of round-off, we could just test whether \texttt{A(i,j) ==0}, 
  % otherwise use \texttt{testzero(A(i,j),norm(A(:,j)),tol)}.
  i = l; while ((i<=n) && testzero(A(i,j),norm(A(:,j)),tol)), i = i+1; end
  if (i<=n)
    % Swap rows \texttt{i} and \texttt{j} in order to move a non-zero matrix entry to
    % position \texttt{(l,j)}. Swapping rows is a valid row transformation!
    A([i,l],:) = A([l,i],:); 
    % Rescale \texttt{l}th row to ensure \texttt{A(l,j) = 1}. This is a row
    % transformation of type (ii)
    A(l,:) = A(l,:)/A(l,j);
    % Now subtract multiples of \texttt{l}th row from the rows below to
    % annihilate the matrix entries in the \texttt{j}th column. This amounts to
    % applying row transformations of type (i).
    for k=[(1:l-1),(l+1:n)],  A(k,:) = A(k,:) - A(k,j)*A(l,:); end
    l = l+1; 
  end
  fprintf('i = %d, l= %d, j = %d\n',i,l,j); A,
  j = j+1;
end
b = A(:,m+1); A = A(:,1:m); 
end

% Surrogate function for testing whether a matrix entry vanishes
function iszero = testzero(z,ref,tol)
  if (nargin == 1), iszero = (z == 0);
  else iszero = (abs(z) <= tol*ref); end
end