function [P,L,U]=plufactor(A)
% PA=LU factorization for n x n A
% L is lower triangular
% U is upper triangular
% WARNING: possible division by zero

n=length(A);
U=A;L=eye(n);P=L;

for p=1:n-1
       [y,k]=max(abs(U(p:n,p)));k=k+p-1;
       aux=U(p,p:n);U(p,p:n)=U(k,p:n);U(k,p:n)=aux;
       aux=L(p,1:p-1);L(p,1:p-1)=L(k,1:p-1);L(k,1:p-1)=aux;
       aux=P(p,:);P(p,:)=P(k,:);P(k,:)=aux;

      for q=p+1:n
              L(q,p)=U(q,p)/U(p,p);
              U(q,p:n)=U(q,p:n)-U(p,p:n)*L(q,p);
      end
end