function [L,U,p]=CroutP(a)
% Factorización LU  de la matriz a por Crout con pivotación
%
 n=size(a,1); p=1:n;
 for k=1:n
    i=k:n; a(i,k)=a(i,k)-a(i,1:k-1)*a(1:k-1,k);
    [r,m]=max(abs(a(k:n,k))); m=m+k-1;
    if a(m,k)==0, continue, end
    if k~=m, a([k m],:)=a([m k],:); p([k m])=p([m k]); end
    i=k+1:n; a(k,i)=(a(k,i)-a(k,1:k-1)*a(1:k-1,i))/a(k,k);
 end
 L=tril(a,0); U=triu(a,1)+eye(n,n);
end
