function [x,res]=SOR(A,b,relax,x0,tol,maxiter)
% SOR R?solution de syst?mes lin?aires par la m?thode de surrelaxation
%   [x,res]=SOR(A,b,relax,x0,tol,maxiter) r?soud Ax=b par la m?thode de
%   SOR avec parametre relax  en utilisant x0 pour choix initial. L'algorithme 
%   s'arr?te quand maxiter op?rations ont ?t? effectu? ou que la norme 2 du 
%   r?sidu devient plus petite que tol. La variable de sortie res est un 
%   vecteur contenant l'historique des normes des r?sidus.

if nargin <3
  error('not enough arguments');
end
if nargin < 6 
  maxiter=500; 
end
if nargin <5
  tol=1e-6;
end
if nargin <4
   x0=zeros(size(A,1),1);
end

L=tril(A,-1);
U=triu(A,1);
D=diag(diag(A));

x=x0;

res(1)=norm(b-A*x);

i=1;
while i<=maxiter & res(i)>tol
  x=(D/relax+L)\(((1-relax)/relax*D-U)*x+b);
  res(i+1)=norm(b-A*x);
  i=i+1;
end

