function varargout = hhess(A) 
%HHESS Householder Hessenberg reduction.
%    [P H] = hhess(A) returns an upper Hessenberg matrix H and
%    orthogonal matrix P such that A = P*H*P'. H and A have
%    the same eigenvalues.
%
%    H = hhess(A) computes the upper Hessenberg matrix orthogonally
%    similar to A. The eigenvalues of H and A are equal.

   n = size(A,1);
   H = A;
	P = eye(n);
   
   for k=1:n-2
      [u, beta] = house(H,k+1,k);
      H(k+1:n,k:n) = H(k+1:n,k:n) - beta*(u*u')*H(k+1:n,k:n);
      P(1:n,k+1:n) = P(1:n,k+1:n) - beta*P(1:n,k+1:n)*(u*u')
      P*A
      H(1:n,k+1:n) = H(1:n,k+1:n) - beta*H(1:n,k+1:n)*(u*u');
      P(1:n,k+1:n) = P(1:n,k+1:n) - beta*P(1:n,k+1:n)*(u*u');
      H(k+2:n,k) = zeros(n-k-1,1);
   end

if nargout == 2
   varargout{1} = P;
   varargout{2} = H;
else
    varargout{1} = H;
end
