function [L, U, P, interchanges] = luhess(H)
%luhess LU factorization of an upper Hessenberg matrix
%using Gaussian elimination with partial pivoting.
%
%   [L U P interchanges] = luhess(H) computes the LU factorization
%   of the upper Hessenberg matrix H. Efficiency is gained because
%   it is only necessary to zero-out one element below each diagonal
%   element.

[m, n] = size(H);
if m ~= n
	error('Matrix H must be square');
end

L = eye(n);
P = eye(n);
interchanges = 0;
	
for i = 1:n-1
   if abs(H(i,i)) == 0 && abs(H(i+1,i)) == 0
      continue;
   end
   if abs(H(i,i)) < abs(H(i+1,i))
      pivotindex = i+1;
		tmp = H(i, i:n);
		H(i, i:n) = H(pivotindex, i:n);
		H(pivotindex, i:n) = tmp;
		tmp = P(i, 1:n);
		P(i, 1:n) = P(pivotindex, 1:n);
		P(pivotindex, 1:n) = tmp;
		tmp = L(i, 1:i-1);
		L(i, 1:i-1) = L(pivotindex, 1:i-1);
		L(pivotindex, 1:i-1) = tmp;
		interchanges = interchanges + 1;
	end
	
	factor = H(i+1,i)/H(i,i);
	H(i+1, i+1:n) = H(i+1, i+1:n) - factor*H(i, i+1:n);
   H(i+1,i)= 0;
	L(i+1,i) = factor;
end

U = triu(H);
