Algorithm for Gauss-Jordan elimination in Matlab

function [A, k] = rref (A, tolerance)

  ## Supress empty list warnings
  eleo = empty_list_elements_ok;
  unwind_protect
    empty_list_elements_ok = 1;
    
    [rows,cols] = size (A);

    if (nargin < 2)
      tolerance = eps * max (rows, cols) * norm (A, inf);
    endif

    used = zeros(1,cols);
    r = 1;
    for c=1:cols
      ## Find the pivot row
      [m, pivot] = max (abs (A (r:rows, c)));
      pivot = r + pivot - 1;
      
      if (m <= tolerance)
      	## Skip column c, making sure the approximately zero terms are
	## actually zero.
      	A (r:rows, c) = zeros (rows-r+1, 1);

      else
	## keep track of bound variables
	used (1, c) = 1;

      	## Swap current row and pivot row
	A ([pivot, r], c:cols) = A ([r, pivot], c:cols);

	## Normalize pivot row
	A (r, c:cols) = A (r, c:cols) / A (r, c);

	## Eliminate the current column
	ridx = [1:r-1, r+1:rows];
	A (ridx, c:cols) = A (ridx, c:cols) - A (ridx, c) * A(r, c:cols);

	## Check if done
      	if (r++ == rows) break; endif
      endif
    endfor
    k = find(used);

  unwind_protect_cleanup
    ## Restore state
    empty_list_elements_ok = eleo;
  end_unwind_protect
endfunction