function [x,iter,relresid] = gausseidel(A,b,x0,tol,maxiter)
%GAUSSEIDEL Gauss-Seidel method for the solution of the
%n x n system Ax = b.
%
%   [x,iter,relresid] = gausseidel(A,b,x0,tol,maxiter) computes
%   the solution x of Ax = b using the Gauss-Seidel iterative method.
%   x0 is the initial approximation to the solution,
%   tol is the error tolerance, and maxiter is the number of
%   iterations to execute. If the Gauss-Seidel method converges, iter
%   contains the number of iterations required to meet the
%   error tolerance. If the Gauss-Seidel iteration does not converge,
%   iter = -1. relresid is the relative residual obtained by
%   the iteration.

[m,n] = size(A);
if m~=n
	disp('matrix A  is not square')  
	return;
end;

k = 1;
x = x0;
for k = 1:maxiter
	x(1) = (1.0/A(1,1))*(b(1) - A(1,2:n)*x(2:n));
	for i = 2:n-1
		x(i) = (1.0/A(i,i))*(b(i) - ... 
				 A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x(i+1:n));
	end
	x(n) = (1.0/A(n,n))*(b(n) - A(n,1:n-1)*x(1:n-1));
   relresid = norm(b-A*x)/norm(b);
	if relresid < tol
		iter = k;
		return;
	end
end

iter = -1;
