function [Q,R,I] = qr_factor(A)

% function [Q,R,I] = qr_factor(A)
% qr-factorization A(:,I) = Q*R 
% with Q unitary and R uppertriangular

[m,n] = size(A); 
A = [A, eye(m)]; I = [1:n];

for i = 1:min(m-1,n);
   J = [i:m]; 
   d = sum(A(J,i:n).^2);
   [d, k] = max(d) ;
   if d == 0; break; end; 
   A(:,[i,k+(i-1)]) = A(:,[k+(i-1),i]); 
   I([i,k+(i-1)]) = I([k+(i-1),i]);
   A(J,i:end)  = h_trf(A(J,i:end));
end;

R = A(:,1:n); Q = A(:,n+1:end)';
