function [x,k,flag] = SOR (A,b,w,x0,toll,kmax)
% SOR Metodo SOR per la risoluzione di un sistema lineare
%     con test di arresto sulla norma della differenza di due 
%     iterate successive
%
% Uso:
%   [x,k,flag] = SOR (A,b,w,x0,toll,kmax)
%
% Dati di ingresso:
%   A    matrice dei coefficienti
%   b    vettore colonna dei termini noti
%   w    fattore di rilassamento
%   x0   vettore colonna iniziale 
%   toll tolleranza per il test di arresto  
%   kmax numero massimo di iterazioni
%
% Dati di uscita:
%   x    array che contiene per colonne le iterate (vettori) del metodo
%   k    numero delle iterazioni effettuate
%   flag vale 0 se e' stato possibile costruire la matrice di iterazione
%        vale 1 se la matrice M risulta singolare

% Inizializzazioni
k = 1;
x(:,1) = x0;
test = toll+1;
flag = 0;

% Controllo sulla singolarita' di M 

% Estrae il vettore corrispondente alla diagonale di A
D = diag(A);

if prod(D) == 0
    disp('*** Matrice D singolare: impossibile costruire la matrice B')
    disp(' ')
    flag = 1;
    return
end

% Costruisce la matrice M = (D-wE)/w

D = diag(D);
M = (D+w*tril(A,-1))/w;

% Costruisce la matrice N = ((1-w)D+wE)/w
N = ((1-w)*D-w*triu(A,1))/w;


% Costruzione matrice di iterazione B e vettore q
B = inv(M)*N; %oppure M\N;
q = inv(M)*b; % oppure M\b;

% Ciclo metodo SOR
while  (test >= toll) && (k < kmax)   
    k = k+1;  
    x(:,k) = B*x(:,k-1)+q;
    test = norm(x(:,k)-x(:,k-1));
end;
k = k-1;