clear
close all

A = mmread('bcsstk01.mtx');
B = mmread('bcsstk06.mtx');
C = mmread('bcsstk14.mtx');
D = mmread('bcsstk15.mtx');
E = mmread('bcsstk17.mtx');

nA = length(A);
nB = length(B);
nC = length(C);
nD = length(D);
nE = length(E);
n  = [nA nB nC nD nE];
nn = [nnz(A) nnz(B) nnz(C) nnz(D) nnz(E)];

%%
xA = ones(nA,1);
xB = ones(nB,1);
xC = ones(nC,1);
xD = ones(nD,1);
xE = ones(nE,1);
bA = A*xA;
bB = B*xB;
bC = C*xC;
bD = D*xD;
bE = E*xE;
temp_pcg   = 0;
temp_bicg  = 0;
temp_gmres = 0;
temp_pcg_p   = 0;
temp_bicg_p  = 0;
temp_gmres_p = 0;

res_pcg      = zeros(length(n),1);
res_bicg     = zeros(length(n),1);
res_gmres    = zeros(length(n),1);
time_pcg     = zeros(length(n),1);
time_bicg    = zeros(length(n),1);
time_gmres   = zeros(length(n),1);
relres_pcg   = zeros(length(n),1);
relres_bicg  = zeros(length(n),1);
relres_gmres = zeros(length(n),1);
iter_pcg     = zeros(length(n),1);
iter_bicg    = zeros(length(n),1);
iter_gmres   = zeros(length(n),1);

res_pcg_p      = zeros(length(n),1);
res_bicg_p     = zeros(length(n),1);
res_gmres_p    = zeros(length(n),1);
time_pcg_p     = zeros(length(n),1);
time_bicg_p    = zeros(length(n),1);
time_gmres_p   = zeros(length(n),1);
relres_pcg_p   = zeros(length(n),1);
relres_bicg_p  = zeros(length(n),1);
relres_gmres_p = zeros(length(n),1);
iter_pcg_p     = zeros(length(n),1);
iter_bicg_p    = zeros(length(n),1);
iter_gmres_p   = zeros(length(n),1);


runs = 5;
for i = 1:runs
    [xA_pcg,flag_A_pcg,relres_A_pcg,iter_A_pcg] = pcg(A,bA,[],10000);
    tic
    [xA_pcg,flag_A_pcg,relres_A_pcg,iter_A_pcg] = pcg(A,bA,[],10000);
    temp = toc;
    temp_pcg = temp_pcg + temp;
    
    [xA_bicg,flag_A_bicg,relres_A_bicg,iter_A_bicg] = bicg(A,bA,[],10000);
    tic
    [xA_bicg,flag_A_bicg,relres_A_bicg,iter_A_bicg] = bicg(A,bA,[],10000);
    temp = toc;
    temp_bicg = temp_bicg + temp;
    
    [xA_gmres,flag_A_gmres,relres_A_gmres,iter_A_gmres] = gmres(A,bA,[],[],size(A,1));
    tic
    [xA_gmres,flag_A_gmres,relres_A_gmres,iter_A_gmres] = gmres(A,bA,[],[],size(A,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
    
    setup.type = 'nofill';
    
    [L,U] = ilu(A,setup);
    [xA_pcg_p,flag_A_pcg_p,relres_A_pcg_p,iter_A_pcg_p] = pcg(A,bA,[],10000,L,U);
    tic
    [L,U] = ilu(A,setup);
    [xA_pcg_p,flag_A_pcg_p,relres_A_pcg_p,iter_A_pcg_p] = pcg(A,bA,[],10000,L,U);
    temp = toc;
    temp_pcg_p = temp_pcg_p + temp;
    
    [L,U] = ilu(A,setup);
    [xA_bicg_p,flag_A_bicg_p,relres_A_bicg_p,iter_A_bicg_p] = bicg(A,bA,[],10000,L,U);
    tic
    [L,U] = ilu(A,setup);
    [xA_bicg_p,flag_A_bicg_p,relres_A_bicg_p,iter_A_bicg_p] = bicg(A,bA,[],10000,L,U);
    temp = toc;
    temp_bicg_p = temp_bicg_p + temp;
    
    [xA_gmres_p,flag_A_gmres_p,relres_A_gmres_p,iter_A_gmres_p] = gmres(A,bA,[],[],size(A,1));
    tic
    [xA_gmres_p,flag_A_gmres_p,relres_A_gmres_p,iter_A_gmres_p] = gmres(A,bA,[],[],size(A,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
end

res_pcg(1)      = norm(xA-xA_pcg);
res_bicg(1)     = norm(xA-xA_bicg);
res_gmres(1)    = norm(xA-xA_gmres);
time_pcg(1)     = temp_pcg/runs;
time_bicg(1)    = temp_bicg/runs;
time_gmres(1)   = temp_gmres/runs;
relres_pcg(1)   = relres_A_pcg;
relres_bicg(1)  = relres_A_bicg;
relres_gmres(1) = relres_A_gmres;
iter_pcg(1)     = iter_A_pcg;
iter_bicg(1)    = iter_A_bicg;
iter_gmres(1)   = max(iter_A_gmres);

res_pcg_p(1)      = norm(xA-xA_pcg_p);
res_bicg_p(1)     = norm(xA-xA_bicg_p);
res_gmres_p(1)    = norm(xA-xA_gmres_p);
time_pcg_p(1)     = temp_pcg_p/runs;
time_bicg_p(1)    = temp_bicg_p/runs;
time_gmres_p(1)   = temp_gmres_p/runs;
relres_pcg_p(1)   = relres_A_pcg_p;
relres_bicg_p(1)  = relres_A_bicg_p;
relres_gmres_p(1) = relres_A_gmres_p;
iter_pcg_p(1)     = iter_A_pcg_p;
iter_bicg_p(1)    = iter_A_bicg_p;
iter_gmres_p(1)   = max(iter_A_gmres_p);
%%

temp_pcg   = 0;
temp_bicg  = 0;
temp_gmres = 0;
temp_pcg_p   = 0;
temp_bicg_p  = 0;
temp_gmres_p = 0;

for i = 1:runs
    [xB_pcg,flag_B_pcg,relres_B_pcg,iter_B_pcg] = pcg(B,bB,[],10000);
    tic
    [xB_pcg,flag_B_pcg,relres_B_pcg,iter_B_pcg] = pcg(B,bB,[],10000);
    temp = toc;
    temp_pcg = temp_pcg + temp;
    
    [xB_bicg,flag_B_bicg,relres_B_bicg,iter_B_bicg] = bicg(B,bB,[],10000);
    tic
    [xB_bicg,flag_B_bicg,relres_B_bicg,iter_B_bicg] = bicg(B,bB,[],10000);
    temp = toc;
    temp_bicg = temp_bicg + temp;
    
    [xB_gmres,flag_B_gmres,relres_B_gmres,iter_B_gmres] = gmres(B,bB,[],[],size(B,1));
    tic
    [xB_gmres,flag_B_gmres,relres_B_gmres,iter_B_gmres] = gmres(B,bB,[],[],size(B,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
    
    setup.type = 'nofill';
    
    [L,U] = ilu(B,setup);
    [xB_pcg_p,flag_B_pcg_p,relres_B_pcg_p,iter_B_pcg_p] = pcg(B,bB,[],10000,L,U);
    tic
    [L,U] = ilu(B,setup);
    [xB_pcg_p,flag_B_pcg_p,relres_B_pcg_p,iter_B_pcg_p] = pcg(B,bB,[],10000,L,U);
    temp = toc;
    temp_pcg_p = temp_pcg_p + temp;
    
    [L,U] = ilu(B,setup);
    [xB_bicg_p,flag_B_bicg_p,relres_B_bicg_p,iter_B_bicg_p] = bicg(B,bB,[],10000,L,U);
    tic
    [L,U] = ilu(B,setup);
    [xB_bicg_p,flag_B_bicg_p,relres_B_bicg_p,iter_B_bicg_p] = bicg(B,bB,[],10000,L,U);
    temp = toc;
    temp_bicg_p = temp_bicg_p + temp;
    
    [xB_gmres_p,flag_B_gmres_p,relres_B_gmres_p,iter_B_gmres_p] = gmres(B,bB,[],[],size(B,1));
    tic
    [xB_gmres_p,flag_B_gmres_p,relres_B_gmres_p,iter_B_gmres_p] = gmres(B,bB,[],[],size(B,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
end

res_pcg(2)      = norm(xB-xB_pcg);
res_bicg(2)     = norm(xB-xB_bicg);
res_gmres(2)    = norm(xB-xB_gmres);
time_pcg(2)     = temp_pcg/runs;
time_bicg(2)    = temp_bicg/runs;
time_gmres(2)   = temp_gmres/runs;
relres_pcg(2)   = relres_B_pcg;
relres_bicg(2)  = relres_B_bicg;
relres_gmres(2) = relres_B_gmres;
iter_pcg(2)     = iter_B_pcg;
iter_bicg(2)    = iter_B_bicg;
iter_gmres(2)   = max(iter_B_gmres);

res_pcg_p(2)      = norm(xB-xB_pcg_p);
res_bicg_p(2)     = norm(xB-xB_bicg_p);
res_gmres_p(2)    = norm(xB-xB_gmres_p);
time_pcg_p(2)     = temp_pcg_p/runs;
time_bicg_p(2)    = temp_bicg_p/runs;
time_gmres_p(2)   = temp_gmres_p/runs;
relres_pcg_p(2)   = relres_B_pcg_p;
relres_bicg_p(2)  = relres_B_bicg_p;
relres_gmres_p(2) = relres_B_gmres_p;
iter_pcg_p(2)     = iter_B_pcg_p;
iter_bicg_p(2)    = iter_B_bicg_p;
iter_gmres_p(2)   = max(iter_B_gmres_p);
%%

temp_pcg   = 0;
temp_bicg  = 0;
temp_gmres = 0;
temp_pcg_p   = 0;
temp_bicg_p  = 0;
temp_gmres_p = 0;

for i = 1:runs
    [xC_pcg,flag_C_pcg,relres_C_pcg,iter_C_pcg] = pcg(C,bC,[],10000);
    tic
    [xC_pcg,flag_C_pcg,relres_C_pcg,iter_C_pcg] = pcg(C,bC,[],10000);
    temp = toc;
    temp_pcg = temp_pcg + temp;
    
    [xC_bicg,flag_C_bicg,relres_C_bicg,iter_C_bicg] = bicg(C,bC,[],10000);
    tic
    [xC_bicg,flag_C_bicg,relres_C_bicg,iter_C_bicg] = bicg(C,bC,[],10000);
    temp = toc;
    temp_bicg = temp_bicg + temp;
    
    [xC_gmres,flag_C_gmres,relres_C_gmres,iter_C_gmres] = gmres(C,bC,[],[],size(C,1));
    tic
    [xC_gmres,flag_C_gmres,relres_C_gmres,iter_C_gmres] = gmres(C,bC,[],[],size(C,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
    
    setup.type = 'nofill';
    
    [L,U] = ilu(C,setup);
    [xC_pcg_p,flag_C_pcg_p,relres_C_pcg_p,iter_C_pcg_p] = pcg(C,bC,[],10000,L,U);
    tic
    [L,U] = ilu(C,setup);
    [xC_pcg_p,flag_C_pcg_p,relres_C_pcg_p,iter_C_pcg_p] = pcg(C,bC,[],10000,L,U);
    temp = toc;
    temp_pcg_p = temp_pcg_p + temp;
    
    [L,U] = ilu(C,setup);
    [xC_bicg_p,flag_C_bicg_p,relres_C_bicg_p,iter_C_bicg_p] = bicg(C,bC,[],10000,L,U);
    tic
    [L,U] = ilu(C,setup);
    [xC_bicg_p,flag_C_bicg_p,relres_C_bicg_p,iter_C_bicg_p] = bicg(C,bC,[],10000,L,U);
    temp = toc;
    temp_bicg_p = temp_bicg_p + temp;
    
    [xC_gmres_p,flag_C_gmres_p,relres_C_gmres_p,iter_C_gmres_p] = gmres(C,bC,[],[],size(C,1));
    tic
    [xC_gmres_p,flag_C_gmres_p,relres_C_gmres_p,iter_C_gmres_p] = gmres(C,bC,[],[],size(C,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
end

res_pcg(3)      = norm(xC-xC_pcg);
res_bicg(3)     = norm(xC-xC_bicg);
res_gmres(3)    = norm(xC-xC_gmres);
time_pcg(3)     = temp_pcg/runs;
time_bicg(3)    = temp_bicg/runs;
time_gmres(3)   = temp_gmres/runs;
relres_pcg(3)   = relres_C_pcg;
relres_bicg(3)  = relres_C_bicg;
relres_gmres(3) = relres_C_gmres;
iter_pcg(3)     = iter_C_pcg;
iter_bicg(3)    = iter_C_bicg;
iter_gmres(3)   = max(iter_C_gmres);

res_pcg_p(3)      = norm(xC-xC_pcg_p);
res_bicg_p(3)     = norm(xC-xC_bicg_p);
res_gmres_p(3)    = norm(xC-xC_gmres_p);
time_pcg_p(3)     = temp_pcg_p/runs;
time_bicg_p(3)    = temp_bicg_p/runs;
time_gmres_p(3)   = temp_gmres_p/runs;
relres_pcg_p(3)   = relres_C_pcg_p;
relres_bicg_p(3)  = relres_C_bicg_p;
relres_gmres_p(3) = relres_C_gmres_p;
iter_pcg_p(3)     = iter_C_pcg_p;
iter_bicg_p(3)    = iter_C_bicg_p;
iter_gmres_p(3)   = max(iter_C_gmres_p);
%%

temp_pcg   = 0;
temp_bicg  = 0;
temp_gmres = 0;
temp_pcg_p   = 0;
temp_bicg_p  = 0;
temp_gmres_p = 0;


for i = 1:runs
    [xD_pcg,flag_D_pcg,relres_D_pcg,iter_D_pcg] = pcg(D,bD,[],10000);
    tic
    [xD_pcg,flag_D_pcg,relres_D_pcg,iter_D_pcg] = pcg(D,bD,[],10000);
    temp = toc;
    temp_pcg = temp_pcg + temp;
    
    [xD_bicg,flag_D_bicg,relres_D_bicg,iter_D_bicg] = bicg(D,bD,[],10000);
    tic
    [xD_bicg,flag_D_bicg,relres_D_bicg,iter_D_bicg] = bicg(D,bD,[],10000);
    temp = toc;
    temp_bicg = temp_bicg + temp;
    
    [xD_gmres,flag_D_gmres,relres_D_gmres,iter_D_gmres] = gmres(D,bD,[],[],size(D,1));
    tic
    [xD_gmres,flag_D_gmres,relres_D_gmres,iter_D_gmres] = gmres(D,bD,[],[],size(D,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
    
    setup.type = 'nofill';
    
    [L,U] = ilu(D,setup);
    [xD_pcg_p,flag_D_pcg_p,relres_D_pcg_p,iter_D_pcg_p] = pcg(D,bD,[],10000,L,U);
    tic
    [L,U] = ilu(D,setup);
    [xD_pcg_p,flag_D_pcg_p,relres_D_pcg_p,iter_D_pcg_p] = pcg(D,bD,[],10000,L,U);
    temp = toc;
    temp_pcg_p = temp_pcg_p + temp;
    
    [L,U] = ilu(D,setup);
    [xD_bicg_p,flag_D_bicg_p,relres_D_bicg_p,iter_D_bicg_p] = bicg(D,bD,[],10000,L,U);
    tic
    [L,U] = ilu(D,setup);
    [xD_bicg_p,flag_D_bicg_p,relres_D_bicg_p,iter_D_bicg_p] = bicg(D,bD,[],10000,L,U);
    temp = toc;
    temp_bicg_p = temp_bicg_p + temp;
    
    [xD_gmres_p,flag_D_gmres_p,relres_D_gmres_p,iter_D_gmres_p] = gmres(D,bD,[],[],size(D,1));
    tic
    [xD_gmres_p,flag_D_gmres_p,relres_D_gmres_p,iter_D_gmres_p] = gmres(D,bD,[],[],size(D,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
end

res_pcg(4)      = norm(xD-xD_pcg);
res_bicg(4)     = norm(xD-xD_bicg);
res_gmres(4)    = norm(xD-xD_gmres);
time_pcg(4)     = temp_pcg/runs;
time_bicg(4)    = temp_bicg/runs;
time_gmres(4)   = temp_gmres/runs;
relres_pcg(4)   = relres_D_pcg;
relres_bicg(4)  = relres_D_bicg;
relres_gmres(4) = relres_D_gmres;
iter_pcg(4)     = iter_D_pcg;
iter_bicg(4)    = iter_D_bicg;
iter_gmres(4)   = max(iter_D_gmres);

res_pcg_p(4)      = norm(xD-xD_pcg_p);
res_bicg_p(4)     = norm(xD-xD_bicg_p);
res_gmres_p(4)    = norm(xD-xD_gmres_p);
time_pcg_p(4)     = temp_pcg_p/runs;
time_bicg_p(4)    = temp_bicg_p/runs;
time_gmres_p(4)   = temp_gmres_p/runs;
relres_pcg_p(4)   = relres_D_pcg_p;
relres_bicg_p(4)  = relres_D_bicg_p;
relres_gmres_p(4) = relres_D_gmres_p;
iter_pcg_p(4)     = iter_D_pcg_p;
iter_bicg_p(4)    = iter_D_bicg_p;
iter_gmres_p(4)   = max(iter_D_gmres_p);

%%
temp_pcg   = 0;
temp_bicg  = 0;
temp_gmres = 0;
temp_pcg_p   = 0;
temp_bicg_p  = 0;
temp_gmres_p = 0;


for i = 1:runs
    [xE_pcg,flag_E_pcg,relres_E_pcg,iter_E_pcg] = pcg(E,bE,[],10000);
    tic
    [xE_pcg,flag_E_pcg,relres_E_pcg,iter_E_pcg] = pcg(E,bE,[],10000);
    temp = toc;
    temp_pcg = temp_pcg + temp;
    
    [xE_bicg,flag_E_bicg,relres_E_bicg,iter_E_bicg] = bicg(E,bE,[],10000);
    tic
    [xE_bicg,flag_E_bicg,relres_E_bicg,iter_E_bicg] = bicg(E,bE,[],10000);
    temp = toc;
    temp_bicg = temp_bicg + temp;
    
    [xE_gmres,flag_E_gmres,relres_E_gmres,iter_E_gmres] = gmres(E,bE,[],[],size(E,1));
    tic
    [xE_gmres,flag_E_gmres,relres_E_gmres,iter_E_gmres] = gmres(E,bE,[],[],size(E,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
    
    setup.type = 'nofill';
    
    [L,U] = ilu(E,setup);
    [xE_pcg_p,flag_E_pcg_p,relres_E_pcg_p,iter_E_pcg_p] = pcg(E,bE,[],10000,L,U);
    tic
    [L,U] = ilu(E,setup);
    [xE_pcg_p,flag_E_pcg_p,relres_E_pcg_p,iter_E_pcg_p] = pcg(E,bE,[],10000,L,U);
    temp = toc;
    temp_pcg_p = temp_pcg_p + temp;
    
    [L,U] = ilu(E,setup);
    [xE_bicg_p,flag_E_bicg_p,relres_E_bicg_p,iter_E_bicg_p] = bicg(E,bE,[],10000,L,U);
    tic
    [L,U] = ilu(E,setup);
    [xE_bicg_p,flag_E_bicg_p,relres_E_bicg_p,iter_E_bicg_p] = bicg(E,bE,[],10000,L,U);
    temp = toc;
    temp_bicg_p = temp_bicg_p + temp;
    
    [xE_gmres_p,flag_E_gmres_p,relres_E_gmres_p,iter_E_gmres_p] = gmres(E,bE,[],[],size(E,1));
    tic
    [xE_gmres_p,flag_E_gmres_p,relres_E_gmres_p,iter_E_gmres_p] = gmres(E,bE,[],[],size(E,1));
    temp = toc;
    temp_gmres = temp_gmres + temp;
end

res_pcg(5)      = norm(xE-xE_pcg);
res_bicg(5)     = norm(xE-xE_bicg);
res_gmres(5)    = norm(xE-xE_gmres);
time_pcg(5)     = temp_pcg/runs;
time_bicg(5)    = temp_bicg/runs;
time_gmres(5)   = temp_gmres/runs;
relres_pcg(5)   = relres_E_pcg;
relres_bicg(5)  = relres_E_bicg;
relres_gmres(5) = relres_E_gmres;
iter_pcg(5)     = iter_E_pcg;
iter_bicg(5)    = iter_E_bicg;
iter_gmres(5)   = max(iter_E_gmres);

res_pcg_p(5)      = norm(xE-xE_pcg_p);
res_bicg_p(5)     = norm(xE-xE_bicg_p);
res_gmres_p(5)    = norm(xE-xE_gmres_p);
time_pcg_p(5)     = temp_pcg_p/runs;
time_bicg_p(5)    = temp_bicg_p/runs;
time_gmres_p(5)   = temp_gmres_p/runs;
relres_pcg_p(5)   = relres_E_pcg_p;
relres_bicg_p(5)  = relres_E_bicg_p;
relres_gmres_p(5) = relres_E_gmres_p;
iter_pcg_p(5)     = iter_E_pcg_p;
iter_bicg_p(5)    = iter_E_bicg_p;
iter_gmres_p(5)   = max(iter_E_gmres_p);

%%
figure(1)
lin1 = loglog(nn,time_pcg,'o-',nn,time_bicg,'o-',nn,time_gmres,'o-');
hold on
lin2 = loglog(nn,time_pcg_p,'x--',nn,time_bicg_p,'x--',nn,time_gmres_p,'x--');
set(lin1,'LineWidth',1.5,'MarkerSize',10);
set(lin2,'LineWidth',1.5,'MarkerSize',10);
set(gca,'FontSize',14,'FontWeight','Bold');
legend('Conjugate Gradient','Biconjugate Gradient','GMRES',...
    'Preconditioned Conj Grad','Preconditioned Biconj','Preconditioned GMRES','Location','NW')
grid on
xlabel('Number of nonzero entries')
ylabel('Time (seconds)')
title('Timings')
hold off
print -dtiff timings

% figure(2)
% lin2 = semilogx(nn,res_pcg./norms','o-',nn,res_bicg./norms','o-',nn,res_gmres./norms','o-');
% set(lin2,'LineWidth',1.5,'MarkerSize',10);
% set(gca,'FontSize',14,'FontWeight','Bold');
% legend('Conjugate Gradient','Biconjugate Gradient','GMRES','Location','NW')
% ylabel('Norm')
% xlabel('Number of nonzero entries')
% title('Compare calculated solution to actual (x-x_{calc})')


figure(4)
clf
lin3 = loglog(nn,iter_pcg,'o-',nn,iter_bicg,'o-',nn,iter_gmres,'o-');
hold on
lin4 = loglog(nn,iter_pcg_p,'x--',nn,iter_bicg_p,'x--',nn,iter_gmres_p,'x--');
set(lin3,'LineWidth',1.5,'MarkerSize',10);
set(lin4,'LineWidth',1.5,'MarkerSize',10);
set(gca,'FontSize',14,'FontWeight','Bold');
legend('Conjugate Gradient','Biconjugate Gradient','GMRES',...
    'Preconditioned Conj Grad','Preconditioned Biconj','Preconditioned GMRES','Location','NW')
grid on
ylabel('Iterations')
xlabel('Number of nonzero entries')
title('Iteration Counts')
hold off
print -dtiff iterations

figure(5)
clf
lin5 = loglog(nn,relres_pcg,'o-',nn,relres_bicg,'o-',nn,relres_gmres,'o-');
hold on
lin6 = loglog(nn,relres_pcg_p,'x--',nn,relres_bicg_p,'x--',nn,relres_gmres_p,'x--');
set(lin5,'LineWidth',1.5,'MarkerSize',10);
set(lin6,'LineWidth',1.5,'MarkerSize',10);
set(gca,'FontSize',14,'FontWeight','Bold');
legend('Conjugate Gradient','Biconjugate Gradient','GMRES',...
    'Preconditioned Conj Grad','Preconditioned Biconj','Preconditioned GMRES','Location','SE')
grid on
ylabel('Residual')
xlabel('Number of nonzero entries')
title('Relative Residuals')
hold off
print -dtiff residuals