function [x, piv, Atilda] = gauss_piv(A,b)
%
%
% [x,piv,Atilda] = Gauss_piv(A,b)
%
% soluzione del sistema lineare Ax=b con il metodo
% di eliminazione di Gauss con pivoting parziale
% A matrice n x n
% b  vettore (termine noto) n x 1
% Atilda  prodotto di L e U
%
% piv rappresenta le eventuali permutazioni
% x vettore soluzione n x 1


pivoting=input('Vuoi effettuare il pivoting parziale? (S/N) ','s') ;
[n,m] = size(A);
if m~=n 
    disp('La matrice non e` quadrata')
    return
end
piv = 1:n;
x=zeros(n,1);
for k=1:n-1
    disp(['PASSO ', num2str(k)])
  % ricerca del pivot
  [~,r] = max(abs(A(k:n,k))) ;
  % r: indice di riga che contiene il pivot
  if pivoting=='S'
  % scambio di righe
    l = r+k-1 ;
    piv([k l]) = piv([l k]) ;
    A([k l],:)=A([l k],:) ;
    b([k l])=b([l k]) ;
  else
    disp('Pivoting non effettuato')        
  end
    pivot=A(k,k) ;
    disp(['pivot=', num2str(pivot)])
  if A(k,k) ~= 0
      A(k+1:n,k) = A(k+1:n,k)/pivot;
      A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - A(k+1:n,k)*A(k,k+1:n);
      b(k+1:n)=b(k+1:n)-A(k+1:n,k)*b(k);
   else
      disp('Matrice singolare!')
      x=[];
      Atilda=[] ;
      return;
  end
end
L=eye(n)+tril(A,-1); % L: contiene i moltiplicatori
disp('matrice L: ')
disp(L)

U=triu(A); %U: risultato dell'eliminazione di Gauss
disp(U)

Atilda=L*U ;
pause


% ==========================
% sostituzione all'indietro
%
x(n)=b(n)/U(n,n) ;
for i=n-1:-1:1
    somma=U(i,i+1:n)*x(i+1:n);
    x(i)=(b(i)-somma)/U(i,i) ;
end
disp('Soluzione: ')
disp(x)
 


% ESEMPI
%
%  A= 1     2     3
%     1     2     4
%    -1     0     6
%  b=[-1 0 6]'
%%
%%
% A= 1    2    3
%    3    2     1
%    2    3     1
% b=[1 0 -1]
%%
%%
% A=[1 1+0.5*10^(-15) 3; 2 2 20; 3 6 4]
