function [x,k] = sor (x,A,b,alpha,tol,m)
%----------------------------------------------------------------
% Description: Solve linear algebraic system, Ax = b, iteratively
%              using the successive relaxation method. 
%
% Inputs:      x     = n by 1 vector containing initial guess
%              A     = n by n coefficient matrix (nonzero diagonal
%                      elements)  
%              b     = n by 1 right-hand side vector
%              alpha = relaxation parameter (alpha > 0)
%              tol   = error tolerance used to terminate search
%              m     = maximum number of iterations (m >= 1)
%
% Outputs:     x = n by 1 solution vector
%              k = number of iterations performed.  If 0 < k < m, 
%                  then the following convergence criterion was 
%                  satisfied where r = b - Ax is the residual error
%                  vector:
%
%                  ||r|| < tol 
%
%                  Zero is returned if one of a diagonal element 
%                  of A was found to be zero.  In this case, the 
%                  equations must be reordered.
%                   
% Notes:       When alpha = 1, the function SR reduces to the Gauss-
%              Seidel method.  If A is a symmetric positive-definite
%              matrix, then the function SR should converge for all
%              relaxation parameters in the following range:
%
%                  0 < alpha < 2 */
%----------------------------------------------------------------

% Initialize
   k = 0;
   n = length(x);
   for i = 1 : n
      if abs(A(i,i)) < eps
          return;
      end
   end
% Iterate
err=0.0;   
for l=1:m
      for i = 1 : n;
         d = A(i,i);
         x(i) = (1 - alpha)*x(i) + alpha*b(i)/d;
         for j = 1 : n
            if j ~= i
               x(i) = x(i) - alpha*A(i,j)*x(j)/d; 
            end
         end
      end
      k = k + 1;
      res=A*x-b
   for i=1:n;
   err=err+res(i)*res(i);
   end
   err=err^0.5
   if err < tol
    disp(k); disp(x)
    return;
   else
   end
   err=0.0;
end
%----------------------------------------------------------------
x=[0 0 0 0]';
A=[-4 1 1 1; 1 -4 1 1; 1 1 -4 1; 1 1 1 -4];
b=[1 1 1 1]';
alpha=1.2;
tol=1.0e-5;
m=1000;
>> sor(x,A,b,alpha,tol,m)