
function [xV_J,xV_GS]=demo_J_GS

example=3;

switch example
    case 1
        A=[11 -5 -5;
            5 12 6;
            6 -4 11];
        
        b=[1 23 13]';
        
        sol=ones(3,1);
        x0=zeros(3,1);
    case 2
        A=gallery('poisson',3);
        sol=ones(size(A,1),1);
        b=A*sol;
        x0=zeros(size(A,1),1);
    case 3
        A=[3 -2 -1; 6 -2 2; -9 7 1];
        b=[0 6 -1]';
        sol=ones(3,1);
        
        % x0=zeros(size(A,1),1);
        x0=(1-10^(-14))*ones(3,1);
end

maxit=200;
tol=10^(-15);

[xV_J, error, iter, flag]  = jacobi(A, x0, b, maxit, tol);

[xV_GS, error, iter, flag]  = sor(A, x0, b, 1, maxit, tol);

eJ=NaN*ones(maxit);
eGS=NaN*ones(maxit);

itJ=size(xV_J,2);

for ii=1:itJ
    eJ(ii)=norm(xV_J(:,ii)-sol,inf);
end

itGS=size(xV_GS,2);

for ii=1:itGS
    eGS(ii)=norm(xV_GS(:,ii)-sol,inf);
end

M=max(itJ,itGS);

for ii=10:20:M
    fprintf('\n \t $%3.0f$ & $%1.1e$ & $%1.1e$ \\\\',ii,eJ(ii),eGS(ii));
end








function [xV, error, iter, flag]  = jacobi(A, x, b, max_it, tol)

%  -- Iterative template routine --
%     Univ. of Tennessee and Oak Ridge National Laboratory
%     October 1, 1993
%     Details of this algorithm are described in "Templates for the
%     Solution of Linear Systems: Building Blocks for Iterative
%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
% [x, error, iter, flag]  = jacobi(A, x, b, max_it, tol)
%
% jacobi.m solves the linear system Ax=b using the Jacobi Method.
%
% input   A        REAL matrix
%         x        REAL initial guess vector
%         b        REAL right hand side vector
%         max_it   INTEGER maximum number of iterations
%         tol      REAL error tolerance
%
% output  x        REAL solution vector
%         error    REAL error norm
%         iter     INTEGER number of iterations performed
%         flag     INTEGER: 0 = solution found to tolerance
%                           1 = no convergence given max_it

iter = 0;                                       % initialization
flag = 0;

bnrm2 = norm( b );
xV=[x];

if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol ) return, end

[m,n]=size(A);
[ M, N ] = split( A , b, 1.0, 1 );              % matrix splitting


for iter = 1:max_it,                            % begin iteration
    
    x_1 = x;
    x   = M \ (N*x + b);                         % update approximation
    xV=[xV x];
    error = norm( x - x_1 ) / norm( x );         % compute error
    if ( error <= tol ), break, end              % check convergence
    
    
    
end

if ( error > tol ) flag = 1; end                % no convergence

% END jacobi.m







function [ M, N, b ] = split( A, b, w, flag )
%
% function [ M, N, b ] = split_matrix( A, b, w, flag )
%
% split.m sets up the matrix splitting for the stationary
% iterative methods: jacobi and sor (gauss-seidel when w = 1.0 )
%
% input   A        DOUBLE PRECISION matrix
%         b        DOUBLE PRECISION right hand side vector (for SOR)
%         w        DOUBLE PRECISION relaxation scalar
%         flag     INTEGER flag for method: 1 = jacobi
%                                           2 = sor
%
% output  M        DOUBLE PRECISION matrix
%         N        DOUBLE PRECISION matrix such that A = M - N
%         b        DOUBLE PRECISION rhs vector ( altered for SOR )

[m,n] = size( A );

if ( flag == 1 ),                   % jacobi splitting
    
    M = diag(diag(A));
    N = diag(diag(A)) - A;
    
elseif ( flag == 2 ),               % sor/gauss-seidel splitting
    
    b = w * b;
    M =  w * tril( A, -1 ) + diag(diag( A ));
    N = -w * triu( A,  1 ) + ( 1.0 - w ) * diag(diag( A ));
    
end;

% END split.m







function [xV, error, iter, flag]  = sor(A, x, b, w, max_it, tol)

%  -- Iterative template routine --
%     Univ. of Tennessee and Oak Ridge National Laboratory
%     October 1, 1993
%     Details of this algorithm are described in "Templates for the
%     Solution of Linear Systems: Building Blocks for Iterative
%     Methods", Barrett, Berry, Chan, Demmel, Donato, Dongarra,
%     Eijkhout, Pozo, Romine, and van der Vorst, SIAM Publications,
%     1993. (ftp netlib2.cs.utk.edu; cd linalg; get templates.ps).
%
% [x, error, iter, flag]  = sor(A, x, b, w, max_it, tol)
%
% sor.m solves the linear system Ax=b using the
% Successive Over-Relaxation Method (Gauss-Seidel method when omega = 1 ).
%
% input   A        REAL matrix
%         x        REAL initial guess vector
%         b        REAL right hand side vector
%         w        REAL relaxation scalar
%         max_it   INTEGER maximum number of iterations
%         tol      REAL error tolerance
%
% output  x        REAL solution vector
%         error    REAL error norm
%         iter     INTEGER number of iterations performed
%         flag     INTEGER: 0 = solution found to tolerance
%                           1 = no convergence given max_it

flag = 0;                                   % initialization
iter = 0;

xV=[x];

bnrm2 = norm( b );
if  ( bnrm2 == 0.0 ), bnrm2 = 1.0; end

r = b - A*x;
error = norm( r ) / bnrm2;
if ( error < tol ) return, end

[ M, N, b ] = split( A, b, w, 2 );          % matrix splitting

for iter = 1:max_it                         % begin iteration
    
    x_1 = x;
    x   = M \ ( N*x + b );                   % update approximation
    xV=[xV x];
    error = norm( x - x_1 ) / norm( x );     % compute error
    if ( error <= tol ), break, end          % check convergence
    
    
end
b = b / w;                                  % restore rhs

if ( error > tol ) flag = 1; end;           % no convergence

% END sor.m
