using LinearAlgebra
A=2*rand(5,5).-1
R=copy(A)
m=size(R)[1]
n=size(R)[2]

H=zeros(m,n-1)
function applyQ(M)
    global n,H
    for j=n-1:-1:1
        v=H[:,j]
        M-=2*v*(v'*M)
    end
    return M
end

for j=1:n-1
    global v,c,R,H
    v=R[:,j]
    v[1:j-1].=0
    c=norm(v)
    if v[j]>0
        v[j]+=c
    else
        v[j]-=c
    end
    v/=norm(v)
    R-=2*v*(v'*R)
    R[j+1:m,j].=0
    H[:,j]=v
end
display(R)
display(H)
println("norm(A-QR)=",norm(A-applyQ(R)))