function [P,L,U]=plufactor(A)
% PA=LU factorization for n x n A
% L is lower triangular
% U is upper triangular

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
