function [x, iter]= itermeth(A,b,x0,nmax,tol,P)
%ITERMETH    Un metodo iterativo generico
%   X = ITERMETH(A,B,X0,NMAX,TOL,P) cerca di risolvere iterativamente
%   il sistema di equazioni lineari A*X=B su X. La matrice A di N-per-N
%   coefficienti deve essere non singolare ed il termine noto B deve
%   avere lunghezza N. Se P='J' viene usato il metodo di Jacobi, se 
%   P='G' viene invece selezionato il metodo di Gauss-Seidel. Altrimenti,
%   P e' una matrice N-per-N non singolare che gioca il ruolo di
%   precondizionatore. Il metodo si arresta quando il rapporto fra la 
%   norma del residuo corrente ed quella del residuo iniziale e' minore
%   di TOL. NMAX prescrive il numero massimo di iterazioni consentite.

%   A.Quarteroni,F.Saleri, Introduzione al Calcolo Scientifico, 2004

[n,n]=size(A); 
if nargin == 6
  if ischar(P)==1
     if P=='J' 
        L = diag(diag(A)); 
        U = eye(n); beta = 1; alpha = 1;
     elseif P == 'G' 
        L = tril(A); 
        U = eye(n); beta = 1; alpha = 1;
     end
  else
     [L,U]=lu(P); beta = 0;
  end
else
  L = eye(n); U = L; beta = 0;
end
iter = 0; 
r = b - A * x0; 
r0 = norm(r); 
err = norm (r); x = x0;
while err > tol & iter < nmax
  iter = iter + 1;    
  z = L\r; 
  z = U\z;   
  if beta == 0
     alpha = z'*r/(z'*A*z);  
  end
  x = x + alpha*z;    
  r = b - A * x;      
  err = norm (r) / r0;
end 