function B=showechelon(A)
%Shows all the steps of Gauss-Jordanian elimination
%with partial row pivoting
%Call: B=showechelon(A)
%Demetrios T. Christopoulos, dchristop@econ.uoa.gr, Spring 2011
format rat
format compact
more off
disp('Let the matrix be:')
A
[m,n] = size(A);
% Compute the numerical tolerance.
tol = max([m,n])*eps*norm(A,'inf');
% Loop over the entire matrix.
i = 1;
j = 1;
k = 0;
while (i <= m) && (j <= n)
   % Find the pivot element of column j.
     [mx,k] = max(abs(A(i:m,j)));
      k = k+i-1;
   if (mx <= tol)
      % The column is negligible, zero it out.
      disp([' The  column ' int2str(j) ' is numerically negligible'])
      A(i:m,j) = zeros(m-i+1,1)
      j = j + 1;
   else
       if (i ~= k)
         % Do partial pivoting: Swap i-th and k-th rows.
         disp([' Do partial row pivoting: Swap rows ' int2str(i) ' and ' int2str(k) ' : ' ])
         disp([' Row(' int2str(i) ') <--> Row(' int2str(k) ') '])
         A([i k],:) = A([k i],:)
      end
     
      %Make the pivot element equal to unity
      %by dividing the pivot row with the pivot element.
      disp([' Make pivot element A(' int2str(i) ',' int2str(j) ') equal to unity:' ])
      disp([' Row(' int2str(i) ') --> Row(' int2str(i) ') / A(' int2str(i) ',' int2str(j) ') '])
      A(i,j:n) = A(i,j:n)/A(i,j)
      %Construct the unit vector at column j 
      %by subtracting proper multiples of the pivot row from all other rows.
      disp([' Do  eliminations in column ' int2str(j) ' : '])
       for k = 1:m
         if k ~= i
            disp([' Row(' int2str(k) ') --> Row(' int2str(k) ') - A(' int2str(k) ',' int2str(j) ') * Row(' int2str(j) ')'])
            A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n)             
         end
      end
      i = i + 1;
      j = j + 1;
   end
end
B=A;