function [x,k] = sr (x,A,b,alpha,tol,m)
%----------------------------------------------------------------
% Usage:       [x,k] = sr (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

   chkvec (x,1,'sr');
   chkvec (b,3,'sr');
   alpha = args (alpha,0,alpha,4,'sr');
   tol   = args (tol,0,tol,5,'sr');   
   m     = args (m,1,m,6,'sr');

   k = 0;
   n = length(x);
   for i = 1 : n
      if abs(A(i,i)) < eps
          return;
      end
   end

% Iterate

   while ((k < m) & (residual(A,b,x) >= tol))
      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;
   end
%----------------------------------------------------------------