function [u, beta] = house(A,i,j)
%HOUSE computes the vector u and the scalar
%beta required to form a Householder reflection that will
%will zero-out A(i+1:n,j).
%
%   [u beta] = house(A,i,j) returns the Householder vector u,
%   and the constant beta required to build the Householder
%   matrix Hu such that Hu*A zeros-out A(i+1:n,j).

   m = size(A,1);

   % find u for the Householder reflection
   u = A(i:m,j); 
   % find the maximum element in absolute value in column j,
   % rows i, ..., m
   colmax = max(abs(u));
   % scale the elements as long as the maximum absolute value
   % is not 0
   if colmax ~= 0
      u = u/colmax;
      colnorm = norm(u);
   else
      colnorm = 0;
   end
   
   if u(1) >= 0
      u(1) = u(1)+ colnorm;
   else
      u(1) = u(1) - colnorm;
   end
   
   if norm(u) ~= 0
      beta = 2/norm(u)^2;
   else
      beta = 0;
   end