frref(A,tol)

function [A,jb] = frref(A,tol)
%FRREF   Fast reduced row echelon form.
%   R = FRREF(A) produces the reduced row echelon form of A.
%   [R,jb] = FRREF(A,TOL) uses the given tolerance in the rank tests.
% 
%   Description: 
%   The algorithm is exactly the same as MATLAB's RREF, only 
%   difference being that vectorization is used instead of for-loop
%   in the original code to subtract multiples of the pivot row from
%   all other rows. A typical speed-up range is about 2-4 times of 
%   the original code. However, the actual speed-up depends on the 
%   size of A. The speed-up is quite considerable (in some cases more
%   than 1000 times) if the number of columns in A far surpass its 
%   number of rows.
%
%   Author: Armin Ataei-Esfahani
%
%   Revisions:
%   25-Sep-2008   Created Function


[m,n] = size(A);

% Compute the default tolerance if none was provided.
if (nargin < 2)
    tol = max(m,n)*eps(class(A))*norm(A,'inf');
    if tol > eps;
        tol = eps(class(A));
    end
end

% Loop over the entire matrix.
i = 1;
j = 1;
jb = [];
% t1 = clock;
while (i <= m) && (j <= n)
   % Find value and index of largest element in the remainder of column j.
   [p,k] = max(abs(A(i:m,j))); k = k+i-1;
   if (p <= tol)
      % The column is negligible, zero it out.
      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else
      % Remember column index
      jb = [jb j];
      % Swap i-th and k-th rows.
      A([i k],j:n) = A([k i],j:n);
      % Divide the pivot row by the pivot element.
      Ai = A(i,j:n)/A(i,j);    
      % Subtract multiples of the pivot row from all the other rows.
      A(:,j:n) = A(:,j:n) - A(:,j)*Ai;
      A(i,j:n) = Ai;
      i = i + 1;
      j = j + 1;
   end
end