function [x,p,L,U] = gauss(A,b)
% Soluzione di un sistema lineare Ax=b
% con metodo di eliminazione di Gauss
% INPUT: A matrice quadrata, b
% OUTPUT: x  soluzione del sistema
%         p  contiene le informazioni sulle eventuali
%            permutazioni
%         A1(=A) matrice uguale al prodotto LU

pivoting=input('Vuoi effettuare il pivoting parziale? (S/N) ','s') ;
[m,n]=size(A) ;
if m~=n
    disp('La matrice non e` quadrata')
    return
end
p=1:n ;
x=zeros(n,1) ;
btilde=b ;

for k=1:n-1
    disp(['STEP ' num2str(k)])
    [~,r]=max(abs(A(k:n,k))) ;  % r: indice che contiene il pivot
    if pivoting=='S'
        l=r+k-1 ;
        p([k l])=p([l k]) ;
        A([k l],:)=A([l k],:) ;
        b([k l])=b([l k]) ;
    else
        disp('Pivoting non effettuato')   
    end
    disp(A)
    % definizione del pivot=A(k,k) ;
    pivot=A(k,k) ; 
    disp(['Elemento pivotale: ' num2str(pivot)])
    disp('Schiaccia un tasto ')
    pause
    if A(k,k)~=0
        A(k+1:n,k)=A(k+1:n,k)/A(k,k) ;
        A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);
        btilde(k+1:n)=btilde(k+1:n)-btilde(k)*A(k+1:n,k) ;
    else
        disp('Elemento pivotale nullo!')
        x=[] ;
        L=[] ;
        U=[] ;
        return
    end
end

L=eye(n)+tril(A,-1) ;
disp('Matrice L:')
disp(L)

U=triu(A) ;
disp(U)
disp(['Termine noto modificato: ' mat2str(btilde)])

% Risoluzione sistema triangolare superiore
% Ux=btilde
x(n)=btilde(n)/U(n,n) ;
for i=n-1:-1:1
    somma=U(i,i+1:n)*x(i+1:n) ;
    x(i)=(btilde(i)-somma)/U(i,i) ;
end
end
