function ShowTwoMethods()
% Compares Jacobi and Gauss-Seidel on random diagonally dominant Ax = b
% problems.

close all

% We will perform this number of steps for each example
nSteps = 100;

% Each example will be a sparse nxn matrix with prescibed density
n = 1000;
density = .01;

for delta = [ 1 .5 .1]
% Set up A with dominance factor = delta,
A = spDD(n,density,delta);

% A = L + D + U where
L = tril(A,-1);
D = diag(diag(A));
U = triu(A,1);

% Jacobi and Gauss-Seidel updates have the form M*xNew = N*xOld + b.
% For Jacobi...
MJ = D;
NJ = -(L+U);

% For Gauss-Seidel...
MG = D+L;
NG = -U;

% Let's set up the rhs b so we know the solution...
xStar = ones(n,1);
b = A*xStar;

% Initializations...
xJ = zeros(n,1);
xG = zeros(n,1);
errJ = zeros(nSteps,1);
errG = zeros(nSteps,1);

% Carry out nSteps for each method...
for k=1:nSteps
xJ = MJ\( NJ*xJ + b);
errJ(k) = norm(xJ-xStar);
xG = MG\( NG*xG + b);
errG(k) = norm(xG-xStar);
end

% The spectral radius = abs(largest eig of M\N) is important..
rhoJ = max(abs(eig(full(MJ\NJ))));
rhoG = max(abs(eig(full(MG\NG))));

% Display how the error is going to zero and the relevance
% of the spectral radius...
figure
semilogy(1:nSteps,errJ,'r',1:nSteps,errG,'b')
ylabel('|| xk - xStar ||','fontsize',16)
xlabel(sprintf('rhoJ = %6.3f rhoG = %6.3f',rhoJ,rhoG),'fontsize',16)
title(sprintf('n = %1d density = %5.3f delta = %5.3f',n,density,delta),'fontsize',16)
set(gca,'fontsize',16)
shg
pause(1)
end

function A = spDD(n,density,delta)
% n is a positive integer, 0 < density < 1, and 0 < delta.
% A is an nxn diagonally dominant sparse matrix.
% Density: nnz(A) is approximately density*n^2
% Diagonal Dominance: abs(A(k,k)) >= delta + abs([A(1:k-1,k);A(k+1:n,k)])
A = sprand(n,n,density);
for k=1:n
A(k,k) = sum(abs(A(:,k)))+delta;
end