%2.3  householder.m

function [U,R] = house(A) % Produce R from Householder reflectors saved in U
[m, n] = size(A); U = zeros(m, n);
for k = 1:n
  w = A(k:m,k);        % start with column k of current A, from diagonal down
  w(1) = w(1)-norm(w); % subtract (w,0,...,0) from a = w. New w = a-r
  u = w/norm(w);       % normalize to unit vector u in the kth reflector H_k
  U(k:m,k) = u;        % save u in U to know the H's that produce Q
  A(k:m,k:n) = A(k:m,k:n)-2*u*(u'*A(k:m,k:n)); % multiply current A by H_k
end
R = triu(A(:,1:n)); % square R from nonzeros on and above diagonal of final A