function fig9_10

%  check CGM for solving Ax=b and
%  try different matrices

n=1000;
fprintf('\n n = %i \n\n',n)

% diagonal matrix
A1=diag(1:n);
cond_A1=cond(A1,2)

% matrix with exactly 3 eigenvalues a, a-1, a+1
a=3;
A2=a*eye(n);
A2(1,n)=1; A2(n,1)=1;
cond_A2=cond(A2,2)

% tri-diagonal
D=2*eye(n,n);
SD1=diag(-ones(n-1,1),-1);
A3=D+SD1+SD1';
cond_A3=cond(A3,2)

% specify solution (xe)
xe=2*(1-0.5*rand(n,1));
%xe=ones(n,1);
b1=A1*xe;
b2=A2*xe;
b3=A3*xe;

% use CGM and return error values
A1_err=cgm(A1,b1,xe);
A2_err=cgm(A2,b2,xe);
A3_err=cgm(A3,b3,xe);

steps1=1:length(A1_err);
steps2=1:length(A2_err);
steps3=1:length(A3_err);
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(steps1,A1_err,'-r','LineWidth',2)
hold on
grid on
set(gca,'MinorGridLineStyle','none')
axis([ 1 1200 1e-12 10])
set(gca,'YTick',[1e-12 1e-8 1e-4 1 ])
hold on
ylabel('Error')
xlabel('Iteration Steps (k)')
%legend(' A_1','Location','South','FontSize',16,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')

loglog(steps2,A2_err,'--k','LineWidth',1.6,'MarkerSize',9)
%legend(' A_1',' A_2','Location','South','FontSize',16,'FontWeight','bold')

loglog(steps3,A3_err,'-b','LineWidth',2.4)
legend(' A_1',' A_2',' A_3','Location','South','FontSize',16,'FontWeight','bold')

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


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

% stopping tolerance
tol=1000*eps;

% starting point
n=length(b);
x=zeros(n,1);

% CGM
r=b-A*x;
d=r;
rr=r'*r;
err=10*tol;
counter=0;
while err>tol && counter<10*n
    counter=counter+1;
    q=A*d;
    alpha=rr/(d'*q);
    bad=alpha*d;
    x=x+bad;
    r=r-alpha*q;
    rr0=rr;
    rr=r'*r;
    beta=rr/rr0;
    d=r+beta*d;
    err=norm(bad,inf);
    error(counter)=norm(x-sol,inf);
    if error(counter)==0
        error(counter)=1e-1000;
    end
end

















