function [ x, k, res ] = preconditioned_gradient( A, b, P, x0, tol, kmax )
% PRECONDITIONED_GRADIENT solve the linear system A x = b by means 
% of the Preconditioned Gadrient method; the preconditioning matrix must be 
% non singular. Stopping criterion based on the residual.
%  [ x, k, res ] = preconditioned_gradient( A, b, P, x0, tol, kmax )
%  Inputs:  A    = matrix (square matrix)
%           b    = vector (right hand side of the linear system)
%           P    = preconditioning matrix (non singular, same size of A)
%           x0   = initial solution (colum vector)
%           tol  = tolerence for the stopping criterion based on residual
%           kmax = maximum number of iterations
%  Outputs: x    = solution vector (column vector)
%           k    = number of iterations at convergence
%           res  = value of the norm of the residual at convergence
%

k = 0;
x = x0; 
r = b - A * x;
res = norm( r );

while( k < kmax && res > tol )        
    z = P \ r;    
    alpha = ( z' * r ) / ( z' * A * z );   
    x = x + alpha * z;
    r = r - alpha * A * z;    
    res = norm( r );
    k = k + 1;        
end

return