function R = ref( A )    % this is a new function for matlab
% REF uses Gaussian elimination with partial pivoting to put A 
%     in row echelon form.  It is called in matlab via the command
%     "ref(X)", or e.g. "ref([1 2 3; 2 3 4; 3 4 5])".
%

[numrows,numcols]=size(A); % Store the size of A

R=A;                     % R will be the result matrix 

row=1;                   % row is the current pivot row
col=1;                   % col is the current pivot column

while (row < numrows) & (col <= numcols)  % loop over rows/cols of matrix

  % First we will do the partial pivoting -- look in the pivot column
  %  for the largest entry and swap that entry's row with the pivot row
  %  this makes sure that our pivots are as big as possible.  We'll also
  %  check to make sure that we found a pivot.

    [maxelement,maxindex] = max( abs( R(row:end,col) ) ); 
                         % maxelement = largest entry in pivot col
                         %      (below the previous pivot)
                         % maxindex = position of maxelement 
                         %      (below the previous pivot)
         
    if (maxelement ~= 0) % if maxelement is not 0, then we found a new pivot

                         % Correct maxindex to be position of maxelement in R
        maxindex = maxindex + row - 1; 
        
                         % Swap the new pivot to be in the correct row
        R([row maxindex],:) = R([maxindex row],:);

                         % Use row operations to kill elements below the pivot
        for kill = (row+1):numrows
            R(kill,col:end) = R(kill,col:end) - R(row,col:end)*R(kill,col)/R(row,col);
        end

  % NOTE: this function would work faster if we replaced the "for" loop above
  %  with a matrix operation.  This is called "vectorizing".
  %  We can vectorize our row reduction by replacing the "for" loop with: 
  %%
  %%    rowOperations = R((row+1):end,col)/R(row,col);  
  %%    R((row+1):end,col:end) = R((row+1):end,col:end) - rowOperations*R(row,col:end);
  %%
  % In the above replacement code, "rowOperations" replaces "R(k,col)/R(row,col)" 
  % row k of "rowOperations*R(row,col:end)" is "R(k,col)/R(row,col)*R(row,col:end)" 
        
        row = row + 1;   % look for the next pivot in the next row

    end                  % end the if statement
    
    col = col + 1;       % look for the next pivot in the next column

end                      % end the while loop

end                      % end the function

