function [L,U,P]=CroutP1(a)
% Factorización L1U de matriz a por Crout con pivotación
%
 n=size(a,1); p=1:n;
 for k=1:n-1
    [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(i,k)=a(i,k)/a(k,k);
    a(i,i)=a(i,i)-a(i,k)*a(k,i);
 end
 L=tril(a,-1)+eye(n,n); U=triu(a); P=eye(n); P=P(p,:);
end
