clear
N = 64;
tol = 1.e-3;

% build linear system
A = delsq(numgrid('S',N+2)); %2D Poisson model problems
n = size(A,1);
b =ones(n,1);

Dv = diag(A); % diagonal part of A as vector
D = diag(Dv); % diagonal part of A as matrix
E = tril(A); % lower triangular part of A

% Jacobi
x = rand(n,1); r = b - A*x;
for i=1:10000,
x = x + r./Dv;
r = b - A*x;
rJ(i)=norm(r)/norm(b);
if rJ(i) < tol, break, end
end
xJ = x;

% Gauss Seidel
x = rand(n,1); r = b - A*x;
for i=1:10000,
x = x + E \ r;
r = b - A*x;
rGS(i) = norm(r)/norm(b);
if rGS(i) < tol, break, end
end
xGS = x;

% SOR
omega = 1+cos(pi/(N+1))/(1+sin(pi/(N+1)))^2; % for SOR
Eomega = ((1-omega)*D + omega*E)/omega; % for SOR

x = rand(n,1); r = b - A*x;
for i=1:10000,
x = x + Eomega \ r;
r = b - A*x;
rSOR(i) = norm(r)/norm(b);
if rSOR(i) < tol, break, end
end
xSOR = x;

%[xCG,flagCG,relresCG,iterCG,rCG] = pcg(A,b,tol,2000);
[xCG,flagCG,relresCG,iterCG] = mycg(A,b,tol,2000);

% Plot relative residuals
figure(1)
semilogy(rJ,'m')
hold on
semilogy(rGS,'r-.')
semilogy(rSOR,'--')
semilogy(relresCG,'r:')
legend('Jacobi','Gauss-Seidel','SOR','myCG')
xlabel('Iterations')
ylabel('Residual norm')
axis tight