function [x,it]=raffinamento(A,b,tol,itmax,pausa)
% RAFFINAMENTO ITERATIVO (applicato alla fattorizzazione LU con pivot in semplice precisione)
% Dati di input
%       A,b: Ax=b sistema che si desidera risolvere
%       tol: precisione richiesta sul residuo b-A*x
%     itmax: numero massimo di iterate consentite
% Dati di output
%         x: soluzione del sistema
%        it: numero di iterazioni effettuate per migliorare la soluzione
%            data dal metodo diretto (fatt. LU in semplice precisione)
%     pausa: 1 se desideriamo mettere in pausa a ogni iterazione, 0
%            altrimenti
% N.B.: A e b dovranno essere passati in doppia precisione
%       in modo da poter calcolare il residuo in doppia precisione
%
% ESEMPIO:
%        n=2000;A=rand(n,n);b=rand(n,1);
%        [x,it]=raffinamento(A,b,1e-10,50,1);
%

if pausa
    pause on
else
    pause off
end
    
As=single(A); % As è in singola precisione
[L,U,P]=lu(As);
% di fatto risolviamo il sistema equivalente P*A*x=P*b
As=P*As;
b=double(P)*b; % double trasforma in doppia precisione una variabile in semplice precisione
A=double(P)*A;
xs=U\(L\b); % soluzione del sistema in semplice precisione
x=double(xs);
r=b-A*x; % residuo in doppia precisione
fprintf('La norma del residuo iniziale è: ||r||=%0.5g \n',norm(r))
disp('Premere un qualsiasi tasto per eseguire un''iterazione del raffinamento')
pause

% raffinamento
it=0;
while norm(r)>=tol && it<itmax
    it=it+1;
    d=U\(L\r); % N.B. d è in semplice precisione
    x=x+double(d); % x è in doppia precisione
    r=b-A*x;
    fprintf('La norma del residuo al passo %d è: ||r||=%0.5g \n',it, norm(r))
    pause
end
