function H = hess1(A)
% HESS1  Reduce a matrix to Hessenberg form.
% (demo only--not efficient)
% Input:
%   A      n-by-n matrix
% Output:
%   H     m-by-n upper Hessenberg form

[n,n] = size(A);
H = A;
for k = 1:n-1
  Q = eye(n);
  z = H(k+1:n,k);
  v = [ -sign(z(1))*norm(z)-z(1); -z(2:end) ];
  P = eye(n-k) - 2*(v*v')/(v'*v);   % HH reflection
  Q(k+1:n,k+1:n) = P;
  H = Q*H*Q';
end
H = triu(H,-1);