function [x,it] = it_method(A,b,type,eeps,maxit,w)
n  = size(A,1);
x  = zeros(n,1);  % x0  then for x^{k+1}
x_ = x;           %          for x^{k}
xgs= x;

if nargin < 5,   w     = 1.5;          end
if nargin < 5,   maxit = 100;          end
if nargin < 4,   eeps  = 1e-4;         end
if nargin < 3,   type  = 'richardson'; end

r0 = b-A*x;   p = r0;
r  = inf*r0;
it = 0;

while (norm(r)/norm(r0) > eeps ) && (it < maxit)
  it = it + 1;
  switch type
    case 'jacobi'
      for i = 1:n,   x(i)   = (1/A(i,i))*( b(i) - A(i,[1:i-1,i+1:end])*x_([1:i-1,i+1:end]) );   end
      x_ = x;
    case 'gauss-seidel'
      for i = 1:n,   x(i)   = (1/A(i,i))*( b(i) - A(i,[1:i-1,i+1:end])*x( [1:i-1,i+1:end]) );   end
    case 'sor'
      for i = 1:n,   xgs(i) = (1/A(i,i))*( b(i) - A(i,[1:i-1,i+1:end])*x( [1:i-1,i+1:end]) );   end
      x  = (1-w)*x_ + w*xgs;
      x_ = x;
    case 'steepest-descent'
      r  = b - A*x;
      al = (r'*r)/(r'*A*r); 
      x  = x + al*r;
    case 'conjugate-gradient'
      r   = b - A*x;
      Ap  = A*p;
      pAp = p'*Ap; 
      al = (r'*p)/pAp;
      x  = x + al*p;
      r  = b - A*x;
      be = -(r'*Ap)/pAp;
      p  = r + be*p;
    otherwise %  'richardson'
      r  = b - A*x;
      x  = x + w*r;
  end
  r  = b - A*x;
end

end