function [A, u] = hzero2(A,i,j,row)
%HZERO2 zeros out all elements A(i+1,j), ..., A(n,j) in
%column j or A(i,j+1), A(i,j+2), ..., A(i,n) in row i of
%an m x n matrix A using a Householder reflection.
%
%   If row = 0, [A u] = hzero2 returns a new m x n matrix A
%   implicitely premultiplied by a Householder matrix Hu.
%   In the return matrix A = Hu*A, all entries in column j
%   below row i are zero.
%
%   If row = 1, [A,u] = hzero2 returns a new m x n matrix A
%   implicitely postmultiplied by a Householder matrix Hu. In
%   the return matrix A = A*Hu, all entries in row i, columns
%   j+1, ..., n are zero.

[m, n] = size(A);

if nargin == 3
   row = 0;
end

% find u for the Householder reflection
if row == 0
   u = A(i:m,j);
else
   u = A(i,j:n);
end

if row == 0
   if norm(u(2:(m-i+1)))< eps
      return;
   end
else
   if norm(u(2:(n-j+1)))< eps
      return;
   end
end

% find the maximum element of u in absolute value
maximum = max(abs(u));
% scale the elements as long as the maximum absolute value
% is not 0
if maximum ~= 0
	u = u/maximum;
	unorm = norm(u);
else
	return;
end
if u(1) >= 0
	u(1) = u(1) + unorm;
else
	u(1) = u(1) - unorm;
end

% implicitely form Hu*A
if norm(u) ~= 0
	beta = 2/norm(u)^2;
else
	beta = 0;
end

if row == 0
   A(i:m,j:n) = A(i:m,j:n) - beta*u*(u'*A(i:m,j:n));
else
   A(i:m,j:n) = A(i:m,j:n) - beta*A(i:m,j:n)*u'*u;
end
   
%    Here is how to build Hu if desired
%    u_tilda = zeros(m,1);
%    u_tilda(i:m) = u;
%    H = eye(m) - beta*u_tilda*u_tilda';
