function [x0,iter] = sucov(A,x0,b,w,ep,numitr);
%SUCOV Successive overrelaxation
%[x,iter] = sucov(A,x0,b,w,ep,numitr) computes the solution x
%of the linear system Ax = b using the successive
%overrelaxation iterative method.  x0 is the initial
%solution, numitr is the number of iterations to be performed,
%specified by the user and w is the relaxation
%parameter. (w > 1). If w = 1 then the successive
%overrelaxation iterative method reduces to the
%Gauss-Seidel iterative method. ep is the tolerance
%If the successive overrelaxation method converged,
%iter contains the iteration number needed to converge. If the successive
%overrelaxation method did not converge, iter contains
%numitr.
%See section 12.2 of the book.
%input  : Matrix A, vectors x0 and b, scalars w and ep and integer numitr
%output : vector x and integer iter

flag = 0;
iter  = 0;
bnorm = norm(b);

if bnorm==0
	bnorm = 1;
end;

r = b-A*x0;
error = norm(r)/bnorm;
 b = w * b;
 M =  w * tril( A, -1 ) + diag(diag( A ));
 N = -w * triu( A,  1 ) + ( 1.0 - w ) * diag(diag( A ));

for iter = 1 : numitr
	x1=x0;
	x0 = M\(N*x0+b);
	error = norm(x0-x1)/norm(x0);
	if ( error <= ep ), break, end          
  end
  b = b / w;