function Y = rref(A)                   

% ----------------------------------------------------------------------
% J. R. Senning
% Gordon College
% February 15, 1995
%
% $Id: rref.m,v 1.4 2005/02/03 18:26:44 senning Exp $
%
% Computes the reduced row echelon form of a rectangular matrix A.
%
% February 10, 1998: Partial pivoting is still used, but now if the
% RREF of the matrix has a row of zeros it is handled properly.  This
% is done by assuming the maxmimal pivot-element candidate is zero if
% it is less than 10 times the machine epsilon.  On a given row, if
% no nonzero pivot element can be found then we give up on that row.
%
% February 3, 2000: Check to make sure that "c" (the pivot column)
% is not greater than the number of columns in the matrix - this can
% happen when the matrix does not have a pivot in each column.
%
% February 3, 2005: Change so the function ends with "endfunction"
# keyword.
% ----------------------------------------------------------------------

my_eps = 10 * eps;

if nargin ~= 1
    disp(setstr(7)), disp('Correct format: Y = rref(A)')
    return
endif

[rows cols] = size(A);

c = 1;
for k = [1 : min(rows, cols)]

    %
    % Find biggest element (in absolute value) in column c between row r
    % and the last row.  If this is zero then we need to move one column
    % to the right and try again.  If we end up at the right-most column
    % and still can't find a nonzero biggest element then we are done.
    %

    r = k;
    biggest = max(abs(A(r:rows,c)));

    while (abs(biggest) < my_eps) && (c < cols)
        c = c + 1;
	biggest = max(abs(A(r:rows,c)));
    endwhile

    %
    % check to see if we are at the end of a row but are still looking
    % for a nonzero pivot.  If we haven't found one then we give up...
    %

    if (abs(biggest) < my_eps) && (c == cols)
	break;
    endif

    %
    % Once a biggest element is found, we need to swap rows so that it
    % will become the pivot element.  First find the row that contains
    % the biggest element then swap the rows.
    %

    while (biggest > abs(A(r,c))) && (r < rows)
	r = r + 1;
    endwhile

    if r != k
	temp_v = A(k,:);
	A(k,:) = A(r,:);
	A(r,:) = temp_v;
    endif

    %
    % Now we can scale the row with the newly-found pivot and use it to
    % eliminate all other entries in the pivot column.
    %

    A(k,:) = A(k,:) / A(k,c);
    for j = [1 : rows]
	if j != k
	    A(j,:) = A(j,:) - A(k,:) * A(j,c);
	endif
    endfor

    %
    % Finally, move to next column (if it exists) and repeat!
    %

    c = c + 1;
    if ( c > cols ) 
	break;
    endif

endfor

%
% Ok, we're done. If the elements in the matrix differ no less than my_eps
% from integer values then we assume that the matrix should have only
% integer values.
%

if (max(abs(A - round(A))) < my_eps)
	Y = round(A);
else
	Y = A;
endif

endfunction