function [H,V]=Hessenberg(A)
% HESSENBERG similarity transformation to Hessenberg form
%   [H,V]=Hessenberg(A) transforms the given matrix A into
%   Hessenberg form H using the orthogonal transformations V,
%   such that H=V'*A*V

n=length(A); V=eye(n); 
for p=2:n-1 
  for q=p+1:n 
    G=planerot(A(p:q-p:q,p-1));                % rotate A(q,p-1)=0
    if G(1,2)~=0                        
      A(p:q-p:q,:)=G*A(p:q-p:q,:);             % A=U^T*A
      A(:,p:q-p:q)=A(:,p:q-p:q)*G';            % A=A*U
      V(:,p:q-p:q)=V(:,p:q-p:q)*G';            % V=V*U
    end
  end
end
H=A;