function [Q,R] = qrHH(A)
% QR Factorization by HH reflectors
% Input  : Q - contains the reflection vectors
% Output : R - R factor of A in the full QR factorization

[n,p] = size(A);

for k = 1:p
    
 %   keyboard
    
    e1 = zeros(n-k+1,1);
    e1(1) = 1;
    qk = (A(k:n,k)  -  norm(A(k:n,k))*e1) / norm(A(k:n,k)  -  norm(A(k:n,k))*e1);
    A(k:n, k:p) =   A(k:n,k:p)  -   (2*qk)*(qk'*A(k:n,k:p));

    Q(1:n-k+1,k) = qk;
end

% keyboard

R = A;

return;

