
%
% Apply the preconditioned CG method to solve the discretized Laplace equation
%
%     - Delta y(x) + react y(x)  = f(x)  x in (0,1)^2,   
%                          y(x)  = 0     x on the boundary
%          
% We use f = 1.
%
% Matthias Heinkenschloss
% March 13, 2009 (last modified March 9, 2017)
%


set(0, 'defaultaxesfontsize',24,'defaultaxeslinewidth',2,...
       'defaultlinelinewidth',2,'defaultpatchlinewidth',2,...
       'defaulttextfontsize',24);



% parameters in the advection diffusion equation
% Laplace
diff  = 1;     % diffusion
adv   = [0;0]; % advection
react = 0;     % reaction

% set mesh points
n1  = 100;
n2  = 100;
h1  = 1/(n1+1);
h2  = 1/(n2+1);
x1  = (0:h1:1)';
x2  = (0:h2:1)';

%% Generate system matrix (lexicographic ordering) using Kronecker products
%  Generate A1
c    = -(diff/h1^2 + adv(1)/h1)*ones(n1,1);
d    = (2*diff/h1^2 + adv(1)/h1 + react)*ones(n1,1);
e    = (-diff/h1^2)*ones(n1,1);
A1   = spdiags([c d e], -1:1, n1, n1);

%  Generate A2
c    = -(diff/h2^2 + adv(2)/h2)*ones(n2,1);
d    = (2*diff/h2^2 + adv(2)/h2 + react)*ones(n2,1);
e    = (-diff/h2^2)*ones(n2,1);
A2   = spdiags([c d e], -1:1, n2, n2);

I1   = speye(n1,n1);
I2   = speye(n2,n2);
A    = kron(I2,A1) + kron(A2,I1);



% right hand side
b = ones(n1*n2,1);


    
% FD  solution
yfd = A\b;

% Get matrics for pointwise preconditioners
E = tril(-A,-1);
F = triu(-A,1);
D = spdiags(diag(A),0,n1*n2,n1*n2);


% solve using CG without preconditioning
[y, res, iter, flag] = mycg( A, b, [] );

figure(1)
fprintf(' CG (no prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res/res(1),'k-'); hold on

% solve using CG with Jacobi preconditioning
[y, res, iter, flag] = mycg( A, b, D );

figure(1)
fprintf(' CG (Jacobi prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res/res(1),'r--'); hold on


% solve using CG with SGS preconditioning
omega = 1;
M     = (1/(omega*(2-omega))) * ((D-omega*E)*(D\(D-omega*F)));

[y, res, iter, flag] = mycg( A, b, M );

figure(1)
fprintf(' CG (SGS prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res/res(1),'b-.'); hold on

% solve using CG with SSOR preconditioning
omega = 1.7;
M     = (1/(omega*(2-omega))) * ((D-omega*E)*(D\(D-omega*F)));


[y, res, iter, flag] = mycg( A, b, M );

figure(1)
fprintf(' CG (SSOR prec.) returned after %d iterations with flag = %d\n', iter,flag)
semilogy((0:iter),res/res(1),'g-'); hold off
xlabel('Iteration'); ylabel('Norm of preconditioned residual')

legend('no prec.', 'Jacobi prec.', 'SGS prec.', 'SSOR prec.',...
       'Location','SouthWest')
axis('tight')

%print -depsc Laplace2D_cg_res.eps

