function [P, L, U] = palu(U)
	n = size(U);
	L = zeros(n);
	P = eye(n);
	for j = 1:n
		%pick max row
		m = j;
		for k = j+1:n
			if abs(U(k,j)) > abs(U(m,j))
				m = k;
			end
		end
		%swap
		v = P(j,:);
		P(j,:) = P(m,:);
		P(m,:) = v;
		v = U(j,:);
		U(j,:) = U(m,:);
		U(m,:) = v;
		U(j+1:n, j) = U(j+1:n, j)./U(j,j);
%		U(j+1:n, j) = zeros(n-j, 1);
		U(j+1:n, j+1:n) = U(j+1:n, j+1:n) - U(j+1:n, j) * U(j, j+1:n);
	end
	for j = 1:n
		L(j,1:j-1) = U(j,1:j-1);
		L(j,j) = 1;
		U(j,1:j-1) = zeros(1,j-1);
	end