function [Q, H] = hess(A)

       [m n] = size(A);
       if m~=n
         disp('square matrix required')
         return
       end
       Q=eye(n); H=A;
       for k=1:n-2
           Y = H(:,k); Z=Y;
           c=sqrt(Y(k+1:n)'*Y(k+1:n));
           Z(k+1)=c;
           Z(k+2:n) = zeros(n-k-1,1);
           P = flip(Y,Z);
           Q = Q*P'; H=P*H*P';
       end