% Jacobi/GS/SOR simple mfile :
% Usage:  x = iter3(A,b,method, TOL,maxit,omega);
%         with method ='jacobi', 'gs', 'sor' // TOL=tolerance for maxit steps
%     or  x = iter3(A,b,'jacobi',1.0E-3,12);    % (no omega) 12 steps
%     or  x = iter3(A,b,'gs', 1.0e-3);          % (no omega) TOL=1.0e-3
%     or  x = iter3(A,b,'gs', 1.0e-3,12,1.25);  % omega=1.25
% e.g.  A = [7     2     0     0
%            2    11    -1     0
%            0    -1     9     3
%            0     0     3    10 ] 
%       b = [ 34 39 18 16 ]' 
%       x = iter3(A,b,'gs')

function x=iter3(A,b,method, TOL,maxit,omega);
if nargin<2, help iter3, return; end
  if nargin<3, method='gs'; end
  if nargin<4, TOL=1.0E-4; end
  if nargin<5, maxit=10;   end
  if nargin<6, omega=1.5;  end
n=size(A,1);
  if n>100, disp('WARNING :  is this too large for a test?'), end
x_exact = A \ b  ;
x=zeros(n,1); x1=ones(n,1); X=[]; k=0;

%============================================================
  L=-tril(A,-1);D=diag(diag(A)); U=-triu(A,1);
 mat_prt(A,D,'|','A','D')
 mat_prt(L,U,':','L','U for A=D - L - U')
if strcmp(method,'jacobi')
disp('_________________________________________ Jacobi')
%%%T_J=D\(L+U), c = D\b  % Optionally
  while k<maxit & norm(x1-x)>TOL | k<1
    x1=x;
    x=D\( (L+U)*x ) + D\b;   X=[X x]; k=k+1;
  end
X
 mat_prt4(X,[],':','All iterates X_Jacobi - matrix version')
 Error = norm(X(:,end)-x_exact) ;

elseif strcmp(method,'gs')
disp('_________________________________________ GS below')
%%T_GS = (D+L)\U, c_gs = (D+L)\b, x=T_GS*x+c_gs; X=[X x]; %% optional

DL=D+L, c=DL\b
  while k<maxit & norm(x1-x)>TOL | k<1
   x1=x;
   x=DL\ (U*x) + c;   X=[X x]; k=k+1;
  end
 mat_prt4(X,[],':','All iterates X_GS - matrix version')
 Error = norm(X(:,end)-x_exact) ;

elseif strcmp(method,'sor') 
disp('_________________________________________ SOR next')
if n<11 & nargin<6
T=D\(L+U);
omega = 2/( 1 + sqrt(1-abs(max(eig(T)))^2 ) ); % optimal omega
end

  [xp kp]=point_sor(A,b,x,maxit,TOL,omega);
  Errorp = norm(xp-x_exact) ;
  fprintf('%s: 2-norm error = %6.1e after %d steps, omega=%f [Pointwise]\n',...
  method, Errorp, kp,omega)

DL=D-omega*L;  %MATRIX version %%%%%%%%%%%%%%%%%%%%
M_SOR =  (1-omega)*D + omega*U ; %% T_SOR=DL\T_SOR 

c_SOR = omega* (DL\b) ;   %<<== Danger of using  omega*DL\b = (omega*DL)\b
  while k<maxit & norm(x1-x)>TOL | k<1
    x1=x;
    x=DL \ (M_SOR*x) + c_SOR; X=[X x]; k=k+1;
  end
  mat_prt4(X,[],':','All iterates X_SOR - matrix version')
  Error = norm(X(:,end)-x_exact) ;
end

fprintf('%s: 2-norm error = %6.1e after %d steps\n',method, Error, k)
  mat_prt4(x,x_exact,'|','x', 'x_exact')
disp('------------------------------------------------------ E N D')
return

function [x,k]=point_sor(A,b,x,maxit,TOL,omega);
disp('_________________________________________ SOR - pointwise')
  C=[];x1=x+1;k=0;
  while k<maxit & norm(x1-x)>TOL | k<1
   x1=x; n=size(A,1);
  for r=1:n % row / eqn no 
    s=0;
  for j=1:(r-1)
    s = s + A(r,j)*x(j);
  end
  for j=(1+r):n
    s = s + A(r,j)*x(j);
  end
    xg = (b(r)-s)/A(r,r);    % Literally GS local
    x(r) = omega* xg + (1-omega)*x(r) ; 
  end 
   C=[C x]; k=k+1; 
 end % k
  mat_prt4(C,[],':','All iterates X_SOR - pointwise')
disp('_________________________________________ SOR - pointwise')
