function [x,itns] = gauss_seidel(A,b,x,tol,maxit)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [x,itns] = gauss_seidel(A,b,x,tol,maxit)
%
% This function uses the Gauss-Seidel method to solve the linear system
% Ax = b. 
%
% Input:
%      A:     Array containing the matrix for the problem.
%      b:     Array containing the right-hand side for the problem.
%      x:     Array containing initial guess for the solution.
%    tol:     The iteration stops with success when the residual
%             norm ||b - A x|| has been reduced by a factor tol.
%  maxit:     Maximum number of iterations.
% Output:
%
%   itns:     Array containing the number of iterations for each method.
%      x:     The approximate solution.
%
% It would be better if the function accessed A in a column-oriented
% fashion rather than row-oriented.
% This implementation is particularly bad since Matlab will
% expand A(i,:) into an n-vector before multiplication.
% See a better implementation in smooth.m, for the multigrid
% project.
%
% Also notice that the termination check doubles the work per iteration.
% It would be better to terminate based on the change in x, but
% this way is better for comparison with the other methods.
%
% Dianne P. O'Leary 04/2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

itns = 0;
r0 = norm(b-A*x);
r = r0;
n = length(x);

while (r/r0 > tol)

   for i=1:n,
     x(i) = (b(i)- A(i,:)*x + A(i,i)*x(i) ) / A(i,i);
   end

   r = norm(b-A*x);
  
   itns = itns + 1;
   if (itns > maxit)
     disp('Gauss-Seidel failed to converge.')
     break
   end
   
end