function [x,L,U,P]=lupivot(A,b)

%------------------------------------
% INPUT:
% A matrice.
% b termine noto.
%
% OUTPUT:
% x soluzione del sistema Ax=b.
% A matrice elaborata dalla funzione.
% r vettore degli scambi.
%-------------------------------------

[n,n]=size(A);
x=zeros(n,1);
y=zeros(n,1);
c=zeros(1,n);
r=1:n;

for p=1:n-1

    % Calcola il pivot relativo alla "p"-sima colonna.
    [max1,j]=max(abs(A(p:n,p)));

    % Scambia righe "p" e "j".
    c=A(p,:);
    A(p,:)=A(j+p-1,:);
    A(j+p-1,:)=c;
    d=r(p);
    r(p)=r(j+p-1);
    r(j+p-1)=d;

    if A(p,p) == 0
       fprintf('\n \t [WARNING]: A NON INVERTIBILE');
       break
    end

    for k=p+1:n
        mult=A(k,p)/A(p,p);
        A(k,p)=mult;
        A(k,p+1:n)=A(k,p+1:n)-mult*A(p,p+1:n);
    end
    
end

% Calcola "y".
y(1)=b(r(1));
for k=2:n
    y(k)=b(r(k))-A(k,1:k-1)*y(1:k-1);
end

% Calcola "x".
x(n)=y(n)/A(n,n);
for k=n-1:-1:1
    x(k)=(y(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end

% Calcoliamo P, L, U.
U=(tril(A'))'; L=A-U+eye(size(A)); 

P=zeros(size(A));
for k=1:n
    P(k,r(k))=1;
end