function [P, Q, L, U] = paqlu(U)
	n = size(U);
	L = zeros(n);
	P = eye(n);
	Q = eye(n);
	for j = 1:n
		%pick max row
		g = j;
		h = j;
		for k = j:n
			for l = j:n
				if abs(U(k,l)) > abs(U(g,h))
					g = k;
					h = l;
				end
			end
		end
		%swap
		v = P(j,:);
		P(j,:) = P(g,:);
		P(g,:) = v;
		v = U(j,:);
		U(j,:) = U(g,:);
		U(g,:) = v;
		v = Q(:,j);
		Q(:,j) = Q(:,h);
		Q(:,h) = v;
		v = U(:,j);
		U(:,j) = U(:,h);
		U(:,h) = 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