% generates Figure 7.4 (left)

% linear system
n = 2000; p = 17389; 
AA = spdiags(primes(p)',0,n,n); e = ones(n,2);
for k = 2.^(0:floor(log2(n))), AA = AA + spdiags(e,[-k k],n,n); end
bb = sparse(n,1); bb(1) = 1;

X = []; N = 10:10:2000;
for n=N
    disp(n);
    A=AA(1:n,1:n); 
    b=bb(1:n);
    % solution by diagonally preconditioned CG
    D=spdiags(diag(A),0,n,n);
    [x,flag,relres,iter] = pcg(A,b,1e-16,100,D);
    X = [X,x(1)];
end

x0=X(end);
err = x0-X';
figure(1);
semilogx(N,-log10(err),'k-','LineWidth',2);
hold on;
for k=4:11
    h=line([2^k,2^k],[1,16]);
    set(h,'LineStyle','-.')
    set(h,'Color','r')
end
semilogx(N,-log10(err),'k-','LineWidth',2);
set(gca,'FontSize',18,'FontWeight','normal','FontName','Times');
xlabel('{\itn}','FontWeight','normal','FontSize',24);
ylabel('Digits','FontWeight','normal','FontSize',24);
axis([10,2000,1,16])
%print -deps2c digitsLog2000.eps