function x=mnhf_GaussSeidel(A,b,x,tol,Nmax)
%MNHF_GAUSSSEIDEL Solve Ax=b using the Gauss-Seidel iterative algorithm.
%
% A - (n-by-n) matrix
% b - (n-by-1) RHS vector
% x - initial guess for the solution, x.
% tol - tolerance
% Nmax - maximum number of iterations
%
% Usage  mnhf_GaussSeidel(A,b,x,1e-8,1000); or
%        mnhf_GaussSeidel(A,b,ones(length(b),1),1e-8,1000);

[r,c] = size(A);       % Compute the size of A.
[rr,cc] = size(b);     % Compute the size of b.
if cc==1
 if isequal(r,c)
  if isequal(c,rr)
   
   count = 0;          % Initialize counter.
   
   while ( (norm(A*x-b)>tol) && (count<=Nmax) )
    % Solve for the elements of x.
    x(1) = 1.0/A(1,1)*(b(1)-A(1,2:r)*x(2:r));
    for ii=2:r-1
     x(ii) = 1.0/A(ii,ii)*(b(ii)-A(ii,1:ii-1)*x(1:ii-1) ...
             -A(ii,ii+1:r)*x(ii+1:r));
    end
    x(r) = 1.0/A(r,r)*(b(r)-A(r,1:r-1)*x(1:r-1));
    count = count+1;  % Increment counter.
   end
   
   if count>Nmax
    error('Maximum number of iterations exceeded.')
   end
   
  else
   error('Size mismatch between matrix and vector.') 
  end
 else
  error('Matrix is not square.')
 end
else
 error('b is not a column vector.')
end

end