function fig9_8

%  laplace's equation and the CGM
%  compare different error measures

% A: from Laplace's equation
N=100;
M=N;
n=N*M;
lambda2=1;
D=2*(1+lambda2)*eye(n,n);
SD1=diag(-lambda2*ones(n-1,1),-1);
for i=N+1:N:n-1
    SD1(i,i-1)=0;
end
SDM=diag(-ones(n-M,1),-M);
A=D+SD1+SD1'+SDM+SDM';
%A=sparse(A);

% specify solution (xe)
xe=2*(1-0.5*rand(n,1));
b=A*xe;

%fprintf('\n n = %i   condition number = %3.1e \n\n',n,cond(A,inf))

% use CGM to solve matrix equation
v=cgm(A,b,xe);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% CGM
function x=cgm(A,b,sol)

% set CGM parameters
tol=1000*eps;

n=length(b);
x=zeros(n,1);
tic
% start iteration
r=b-A*x;
d=r;
rr=r'*r;
err=10*tol;
counter=0;
while err>tol
    counter=counter+1;
    q=A*d;
    alpha=rr/(d'*q);
    x=x+alpha*d;
    r=r-alpha*q;
    rr0=rr;
    rr=r'*r;
    beta=rr/rr0;
    d=r+beta*d;
    iter(counter)=counter;
    error2(counter)=norm(alpha*d,inf);
    error3(counter)=norm(r,inf);
    error(counter)=norm(x-sol,inf);
    err=error(counter);
end
toc

fprintf('\n  Number of CGM Iterations = %i     Error =  %e    Time = %e \n\n',counter,error2(counter),toc)

clf
% get(gcf)
set(gcf,'Position', [25 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.18,'MR',-0.02,'ML',0.06,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)

loglog(iter,error,'--k','LineWidth',1.6)
hold on
loglog(iter,error2,'-r','LineWidth',1.6)
loglog(iter,error3,'-b','LineWidth',1.6)

grid on
set(gca,'MinorGridLineStyle','none')
axis([ 1 1000 1e-12 10])
set(gca,'YTick',[1e-12 1e-8 1e-4 1 ])
ylabel('Error')
xlabel('Iteration Steps')
legend({' Error',' Iteration Error',' Residual'},'Location','SouthWest','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')
hold off

%exportgraphics(gcf,'/Users/mark/Desktop/cgm2.eps')













