clear
N = 128;
tol = 1.0e-6;

% build linear system of Example 2.31 (page 133)
%A = sparse(diag(sqrt(1:N))+diag(cos(1:(N-10)),10)+diag(cos(1:(N-10)),-10));
A = delsq(numgrid('S',N+2)); %2D Poisson model problems
n = size(A,1);
b = A*ones(n,1);

Afun=@(x) A*x;%anonymous function to compute A*x
maxit=1000;
%GMRES without preconditioner
[xCG,flagCG,~,iterCG,relresCG] = gmres(Afun,b,[],tol,maxit);

D=spdiags(diag(A),0,n,n); L=tril(A,-1); U=triu(A,1);
Mfun=@(r) D\r; %Jacobi preconditioner
[xPCG,flagPCG,~,iterPCG,relresPCG] = gmres(Afun,b,[],tol,maxit,Mfun);

omega=1;
%SSOR preconditioner (omega=1 is Gauss-Seidel)
Mfun=@(r) (D+omega*U)\(D*((D+omega*L)\r));
[xPCG1,flagPCG1,~,iterPCG1,relresPCG1] = gmres(Afun,b,[],tol,maxit,Mfun);

omega = 2/(1+sin(pi/(N+1))); % optimal omega for SOR
%Symmetric SOR preconditioner
Mfun=@(r) (D+omega*U)\(D*((D+omega*L)\r));
[xPCG2,flagPCG2,~,iterPCG2,relresPCG2] = gmres(Afun,b,[],tol,maxit,Mfun);

%Incomplete Cholesky Factorization
R = ichol(A); %A approximately equal to R*R'
%norm(A-(R*R').*spones(A),'fro')
Mfun=@(r) (R')\(R\r);
[xPCG3,flagPCG3,~,iterPCG3,relresPCG3] = gmres(Afun,b,[],tol,maxit,Mfun);

% Plot relative residuals
figure(1)
semilogy(relresCG,'r-');hold on
semilogy(relresPCG,'r:');hold on
semilogy(relresPCG1,'r-.');hold on
semilogy(relresPCG2,'r--');hold on
semilogy(relresPCG3,'r-');hold on

legend('GMRES','Pre-GMRES: Jacobi','Pre-GMRES: Gauss-Seidel','Pre-GMRES: SSOR(optimal {\omega})','Pre-GMRES: Incomplete Cholesky')
xlabel('Iterations')
ylabel('Residual norm')