  function [alpha,piv,U] = HessLU(A)
% A is an n-by-n upper Hessenberg matrix.
% alpha and piv are 1-by(n-1) with the property that
%         M_{n-1}*P_{n-1}*...*M_{1}*P_{1}*A = U
% is upper triangular where 
%      P_{k} is the identity if piv(k) = 0 and the identity with rows
%         k and k+1 swapped if piv(k)=1.
%      M_{k} is the identity everywhere except that its (k+1,k) entry 
%         is -alpha(k).

  [n,n] = size(A); alpha = zeros(1,n-1); piv = zeros(1,n-1);
  for k=1:n-1
     % Compute the sole multiplier...
     if abs(A(k+1,k))>abs(A(k,k))
        % Interchange rows k and k+1...
        piv(k) = 1;
        A([k k+1],k:n) = A([k+1 k],k:n);
     end
     alpha(k) = A(k+1,k)/A(k,k);
     % zero A(k+1,k)...
     A(k+1,k:n) = A(k+1,k:n) - alpha(k)*A(k,k:n);
  end
  U = triu(A);