function B = rowred(A)
% rowred(A) returns the row reduced echolon form of a matrix A

[m,n]=size(A);  % A is an (m x n)-matrix
r=1;            % r is used to remember the row of A
k=1;


for i=1:n       % i runs through the columns of A                              

     k = r;     
     while (k<m+1)  & (abs(A(k,i))<0.00001)   
      k = k+1;                                
     end                                      

% k is the row of the first non-zero entry in the i-th column 
% below the r-th row. if there is no such entry, k equals m+1
% if k is smaller then m+1, exchange k-th and r-th row and kill all entries 
% A(r+1,i) up to A(m,i)
     
     if k<m+1                                 
         A([r,k],:)=A([k,r],:);
	     for j=r+1:m 
		     A(j,:) = A(j,:)-A(j,i)/A(r,i)*A(r,:);
	      end
 	  r=r+1;
      end
end

% now A is in echelon form
% below, we reduce A to row-reduced echelon form 


for i=m:-1:1                                 % i runs through the rows

     k = 1; 

     while (k<n+1)  & (abs(A(i,k))<0.00001)  % look for pivot in the i-th row
        k = k+1; 
     end

% if there is a pivot in the i-th row, k is the column of that pivot. 
% otherwise k equals n+1. The loop below kills all the other entries in the 
% k-th column and normalizes the pivot

     if k<n+1
     A(i,:)=1/A(i,k)*A(i,:);
         for j = i-1:-1:1
            A(j,:) = A(j,:)-A(j,k)*A(i,:);
         end
     end
end

%now A is in row-reduced echelon form

B=A;