function [A,b,r,I] = GE_theoretical(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$.
% Return values: A = transformed matrix, b = transformed right hand side,
% r = rank of matrix A, I = indices of columns of transformed matrix agreeing with unit vectors 

[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, so that this
% implementation qualifies as a \textbf{theoretical algorithm} only. 
% It \textbf{must not} be used for numerical computations. 

% First trick: add $\Vb$ as rightmost column to $\VA$ in order to facilitate
% implementation of simultaneous row transformations of $\VA$ and $\Vb$
fprintf('Original %d\\times%d-matrix:\n',n,m+1); A = [A,b], 

% \texttt{j} is the current column index, \texttt{l} the current row index.
j = 1; l = 1; I = [];
while (j<=m)
  % Search from top of column \texttt{j} to find first non-zero entry.
  i = l; while ((i<=n) && A(i,j) == 0), 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!
    fprintf('Swapping rows %d and %d:\n',i,l);
    A([i,l],:) = A([l,i],:); A,
    % Rescale \texttt{l}th row to ensure \texttt{A(l,j) = 1}. This is a row
    % transformation of type (ii)
    fprintf('Recaling row %d:\n',l);
    A(l,:) = A(l,:)/A(l,j); A, I = [I,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)]
      fprintf('Subtracting %f*row %d from row %d:\n',A(k,j),l,k);
      A(k,:) = A(k,:) - A(k,j)*A(l,:); A,
    end
    l = l+1; 
  else
    fprintf('No nonzero entry in column %d below row %d\n',j,l);
  end
  j = j+1;
end
r = l-1;
fprintf('Final reduced step form of the matrix, rank = %d\n',r); A,
b = A(:,m+1); A = A(:,1:m); 
end
