GaussSeidelSOR(A,b,omega,es,maxit)

function [x,ea,iter] = GaussSeidelSOR(A,b,omega,es,maxit)
% GaussSeidel: Gauss Seidel method
%
% x = GaussSeidel(A,b): Gauss Seidel without relaxation
%
% x = GaussSeidel(A,b,omega): Gauss Seidel with relaxation
%
% INPUTS:
%           A = coefficient matrix
%           b = right hand side vector
%           omega = weighting factor
%           es = stop criterion (default = 0.00001%)
%           maxit = max iterations (default = 50)
%
% OUTPUT:
%           x = solution vector
%           ea = approximate error
%           iter = iteration completed

% Implemented by ASHISH MESHRAM
% [e-mail],http://www.facebook.com/ashishmeet

if nargin<2
    error('at least 2 input arguments required');
end
if nargin<5||isempty(maxit)
    maxit=50;
end
if nargin<4||isempty(es)
    es=0.00001;
end
if nargin<3||isempty(omega)
    omega=1;
end

[m,n] = size(A); %---Evaluate the dimesnion of a matrix
 %---Checking input matrix whether its a square matrix or not
if m~=n, error('Matrix A must be square'); end

C = A;

for i = 1:n
    C(i,i) = 0; %---Placing zero on principal diagonal
    x(i) = 0; %---Initializing initial guess by setting it to be 0
end

x = x'; %---Transposing

for i = 1:n
    C(i,1:n) = C(i,1:n)/A(i,i);
end

for i = 1:n
    d(i) = b(i)/A(i,i);
end

iter = 0; %---Initializing counter variable "iter"
while iter<=maxit %---condition for while till counter variable 
    %---is less than  or  equal  to maximum  iteration
    xold = x;
    for i = 1:n
        x(i) = d(i)-C(i,:)*x;
        x(i) = omega*x(i) + (1-omega)*xold(i);
        %---Checking convergence by calculating approximate error
        if x(i) ~= 0
            ea(i) = abs((x(i) - xold(i))/x(i)) * 100;
        end
    end
    iter = iter+1; %---incrementing counter variable
    
    %---exiting while loop if maximum of ea is less than or equal to es
      if max(ea)<=es
        break
      end
%     iter
%     x  
%     ea
end