x = [1; 2; 3; 4]
s = norm(x)
z = [-s; 0; 0; 0]
u = x-z;
u = u/norm(u)
x - 2*u*u'*x;
%
%  test housevec.m for a QR decomposition
%
A = rand(7,3)            % the matrix to be decomposed
u1 = housevec(A(:,1))    % the first Householder vector, from the 1st col of A
A1 = A - 2*u1*u1'*A      % the 1st Householder tranformation from left
H = eye(7)-2*u1*u1';     % the 1st Hous. transf. on I from right
u2 = housevec(A1(2:7,2));   % the 2nd Hous. vec.
A2 = A1;
A2(2:7,:) = A1(2:7,:) - 2*u2*u2'*A1(2:7,:)  % the 2nd H.T. on A1(2:7,:)
H(:,2:7) = H(:,2:7)-H(:,2:7)*2*u2*u2';      % 2nd H.T. on H1(:,2:7)
u3 = housevec(A2(3:7,3));
A3 = A2;
A3(3:7,:) = A2(3:7,:)-2*u3*u3'*A2(3:7,:)    % 3rd H.T. on A2(3:7,:)
H(:,3:7) = H(:,3:7) - H(:,3:7)*2*u3*u3';    % 3rd H.T. on H2(:,3:7)
Q = H
R = A3