function [A] = myrref(A, tolerance)
%myrref.m   Reduced Row Echelon Form
%
%   R = myrref(A) produces the reduced row echelon form of matrix A
%   tolerance default value: max(size(A)) * eps * infinity-norm(A)
%   eps is distance from 1 to next larger floating-point number: 2^{-52} 
%
%  RREF algorithm (basic elements only)
%
% Linear Algebra: Theory and Applications
% Ward Cheney & David Kincaid 
% JBPub.com (c) 2008
%
% Based on notes of Carl de Boor 
% GNU General Public License
%
[m, n] = size(A);
%
% If none provided, compute default tolerance 
%
% number of function input arguments: nargin
%
if (nargin < 2)
   tolerance = max(m, n)*eps*norm(A, 'inf'); 
end%if
%
% Loop over entire matrix A
%
i = 1;
j = 1;
%
%
while  (i <= m)  &  (j <= n)
   %
   % In remainder of column j,  find value and index of largest element 
   %
   [p, k] = max(abs(A(i:m, j)));  
   k = k + i - 1;
   %
   if (p <= tolerance)
      %
      % If column is negligible,  zero it out
      %
      A(i:m, j) = zeros(m - i + 1, 1);
      j = j + 1;
      %
   else
      %
      % Swap rows:  i-th  and  k-th 
      %
      A([i k], j:n) = A([k  i], j:n);
      %
      % Divide pivot row by pivot element
      %
      A(i, j:n) = A(i, j:n)/A(i, j);
      %
      % Subtract multiples of pivot row from all other rows
      %
      for k = [1:i-1  i+1:m]
         A(k, j:n) = A(k, j:n) - A(k, j)*A(i, j:n);
      end%for
      %
      i = i + 1;
      j = j + 1;
      %   
   end%if
   %
end%while
%
end%function