function [x,iter,res] = gaussseidel(A,b,x0,tol,itmax)
% Usage: [x,iter,res] = gaussseidel(A,b,x0,tol,itmax)
% Uses Gauss Seidel iteration to solve Ax = b
% Input variables:
% A = matrix
% b = rhs
% x0 = initial geuss
% tol = l2 tolerance on residue
% itmax = maximum allowable number of iterations
% Output variables:
% x = solution
% iter = number of iterations
% res = l2 residual error
%
% Author: Prof. J. Alam
% Course: AMAT3132 - Numerical Analysis I
% Dept of Math. and Stat., Memorial University.


[n,m] = size(A); % Find size of matrix
if abs(n-m)~0;
    fprintf('ERROR: input A must be a square matrix\n%s');
    return;
end

x = x0; % Initial guess
k = 0; % Initialize iteration counter
error = norm(A*x-b,2); % Initial error

while error > tol & k+1 <= itmax
      k = k+1;
      for i = 1:n
	x(i) = (b(i) - A(i,1:i-1)*x(1:i-1) - ... % Use new values of x
	  A(i,i+1:n)*x(i+1:n))/A(i,i);
      end
      error = norm(A*x-b,2);
      res(k)=error;
end
if(k < itmax);
    fprintf('Gauss-Seidel method has been converged\n%s');
else
    fprintf('Tolerance is NOT satisfied\n%s');
end
iter = k;
res=res(1:iter);