% testGM1p: GMRES with pre-conditioner

n = input('n = ');

for s = 1:4 % s-loop
fprintf('\n\t ------- s = %g -------\n',s);

A = spdiags((n-1:-1:1)',-1,n,n);
A(:,n) = (n:-1:1)'; e = ones(n,1);
D = spdiags(e,s+2,n,n);
A = A + D + D';
b = zeros(n,1); b(1) = 1;

% pre-conditioners
[L,U] = luinc(A,1);
       
% plot eigenvalues
eval(['subplot(22' int2str(s) ');']);
eA = eig(full(A)); 
plot(real(eA),imag(eA),'b.'); hold on
eB = eig(full(U\(L\A))); 
plot(real(eB),imag(eB),'r*'); hold off;
axis('square'); drawnow;

bnrm = norm(b); tol = 1e-12;
if exist('myGMRES1p');
    fprintf('myGMRES1p:\n');
    t0 = cputime;
    [x,iter] = myGMRES1p(A,b,tol);
    cput = cputime - t0;
    rres = norm(A*x - b)/(1 + bnrm);
    fprintf('          Without pre-conditioner\n');
    fprintf(' Iter %3i: rel_res = %6.2e  cpu = %6.2e\n',...
        iter,rres,cput);
    t0 = cputime;
    [x,iter] = myGMRES1p(A,b,tol,L,U);
    cput = cputime - t0;
    rres = norm(A*x - b)/(1 + bnrm);
    fprintf('          With pre-conditioner\n');
    fprintf(' Iter %3i: rel_res = %6.2e  cpu = %6.2e\n',...
        iter,rres,cput);
end

if exist('zGMRES1p');
    fprintf(' zGMRES1p:\n');
    t0 = cputime;
    [x,iter] = zGMRES1p(A,b,tol);
    cput = cputime - t0;
    rres = norm(A*x - b)/(1 + bnrm);
    fprintf('          Without pre-conditioner\n');
    fprintf(' Iter %3i: rel_res = %6.2e  cpu = %6.2e\n',...
        iter,rres,cput);
    t0 = cputime;
    [x,iter] = zGMRES1p(A,b,tol,L,U);
    cput = cputime - t0;
    rres = norm(A*x - b)/(1 + bnrm);
    fprintf('          With pre-conditioner\n');
    fprintf(' Iter %3i: rel_res = %6.2e  cpu = %6.2e\n',...
        iter,rres,cput);
end

end % s-loop