function [Q,R]=qr_house2(A)
% QR_HOUSE Householder QR decomposition.
% Given an (m,n) matrix A. This function computes the QR decomposition A=QR
% using Algorithm 5.2.1 in Golub and Van Loan
% [Q,R]=qr_house2(A)
% On input
% A - An (m x n) matrix with m >=n
% On output
% R - The (m x n) upper triangular matrix R in QR decomposition, QR=A
% Q - The orthogonal matrix in the QR decomposition such that Q'*Q=I

% Author: David Di Ruscio

[m,n]=size(A);

Q=eye(m,m);
% Caclculate upper triangular matrix R and A overwritten by this matrix.
for j=1:n
    [v,beta]=house(A(j:m,j));
    
    % Recover orthogonal Q
    H=(eye(m-j+1)-beta*v*v'); 
    Ht=eye(m,m); Ht(j:m,j:m)=H;
    Q=Ht*Q;
    % End recover Q
    
    %A(j:m,j:n)=(eye(m-j+1)-beta*v*v')*A(j:m,j:n)    
    A(j:m,j:n)=H*A(j:m,j:n);
    if j < m
        A(j+1:m,j)=v(2:m-j+1);
    end
    
end
Q=Q';
R=triu(A); 

