function fig9_12

%  comparing PCGM and CGM:  figures 9.12 and 9.13

icase=2

if icase==2
    % 100x100 random
    n=100
    [Q,R]=qr(randn(n,n)); D=diag(logspace(0,4,n));
    A=Q'*D*Q;
    Imax=2*n;
elseif icase==1
    % tri-diagonal
    n=1000;
    D=2*eye(n,n);
    SD1=diag(-ones(n-1,1),-1);
    A=D+SD1+SD1';
    %A=sparse(A);
    Imax=n+20;
end

cond_A=cond(A,inf)
sol=rand(n,1);
sol=2*(1-0.5*rand(n,1));
b=A*sol;

% use SSOR PCGM
LM=tril(A)+diag(diag(A));
UM=eye(n)+inv(diag(diag(A)))*tril(A)';
M=LM*UM;
tic
%[iters,errors,errorIs]=pcgm(A,b,M,sol);
[iters,errors,errorIs]=pcgmLU(A,b,LM,UM,sol);
ssor_time=toc
B=chol(inv(M));
cond_BAB=cond(B*A*B',2)

% use CGM to solve matrix equation
tic
[iter,error,errorI]=cgm(A,b,sol);
cgm_time=toc

figure(1)
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)

semilogy(iter,error,'-b','LineWidth',1.6)
hold on

grid on
set(gca,'MinorGridLineStyle','none')
axis([ 1 Imax 1e-12 4])
set(gca,'YTick',[1e-12 1e-8 1e-4 1 ])
ylabel('Error')
xlabel('Iteration Steps')
%legend(' CGM','Location','SouthWest','FontSize',14,'FontWeight','bold');
set(gca,'FontSize',16,'FontWeight','bold')

figure(2)
clf
% get(gcf)
set(gcf,'Position', [19 794 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)

semilogy(iter,errorI,'-b','LineWidth',1.6)
hold on

grid on
set(gca,'MinorGridLineStyle','none')
axis([ 1 Imax 1e-12 4])
set(gca,'YTick',[1e-12 1e-8 1e-4 1 ])
ylabel('Iterative Error')
xlabel('Iteration Steps')
%legend(' CGM','Location','SouthWest','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',16,'FontWeight','bold')

%pause

figure(1)
semilogy(iters,errors,'-r','LineWidth',1.6)
legend({' CGM',' PCGM'},'Location','SouthWest','FontSize',16,'FontWeight','bold')

% if icase==1
%     exportgraphics(gcf,'/Users/mark/Desktop/pcgmA.eps')
% elseif icase==2
%     exportgraphics(gcf,'/Users/mark/Desktop/pcgm1.eps')
% end

figure(2)
semilogy(iters,errorIs,'-r','LineWidth',1.6)
legend(' CGM',' PCGM','Location','SouthWest','FontSize',16,'FontWeight','bold')

% if icase==1
%     exportgraphics(gcf,'/Users/mark/Desktop/pcgmB.eps')
% elseif icase==2
%     exportgraphics(gcf,'/Users/mark/Desktop/pcgm2.eps')
% end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction for the PCGM
% uses M (not factor of M)
function [iter,error,errorI]=pcgm(A,b,M,sol)

% set PCGM parameters
tol=1e-12;

n=length(b);
x=zeros(n,1);

M=sparse(M);

% start iteration
r=b-A*x;
z=M\r;
d=z;
rz=r'*z;
err=10*tol;
counter=0;
while err>tol && counter<2*n
    counter=counter+1;
    q=A*d;
    alpha=rz/(d'*q);
    x=x+alpha*d;
    r=r-alpha*q;
    z=M\r;
    rz0=rz;
    rz=r'*z;
    beta=rz/rz0;
    d=z+beta*d;
    iter(counter)=counter;
    error(counter)=norm(x-sol,inf);
    errorI(counter)=norm(alpha*d,inf);
    err=error(counter);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction for the PCGM
% uses M (not factor of M)
function [iter,error,errorI]=pcgmLU(A,b,L,U,sol)

% set PCGM parameters
tol=1e-12;

n=length(b);
x=zeros(n,1);

% start iteration
r=b-A*x;
%z=M\r;
y=L\r;
z=U\y;
d=z;
rz=r'*z;
err=10*tol;
counter=0;
while err>tol && counter<2*n
    counter=counter+1;
    q=A*d;
    alpha=rz/(d'*q);
    x=x+alpha*d;
    r=r-alpha*q;
    %z=M\r;
    y=L\r;
    z=U\y;
    rz0=rz;
    rz=r'*z;
    beta=rz/rz0;
    d=z+beta*d;
    iter(counter)=counter;
    error(counter)=norm(x-sol,inf);
    errorI(counter)=norm(alpha*d,inf);
    err=error(counter);
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% subfunction for the CGM
function [iter,error,errorI]=cgm(A,b,sol)

% set CGM parameters
tol=1e-12;

n=length(b);
x=zeros(n,1);

% start iteration
r=b-A*x;
d=r;
rr=r'*r;
err=10*tol;
counter=0;
while err>tol && counter<2*n
    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;
    error(counter)=norm(x-sol,inf);
    errorI(counter)=norm(alpha*d,inf);
    err=error(counter);
end

















