% generates Figure 7.4 (right)

% linear system
n = 20000; p = 224737; 
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 = 100:100:20000;
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'; err = max(err,eps);
figure(1);
semilogx(N,-log10(err),'k-','Linewidth',2);
hold on
set(gca,'FontSize',18,'FontWeight','normal','FontName','Times');
xlabel('{\itn}','FontWeight','normal','FontSize',24);
ylabel('Digits','FontWeight','normal','FontSize',24);
axis([800,21000,2,16])
hold on
for k=7:14
    h=line([2^k,2^k],[2,16]);
    set(h,'LineStyle','-.')
    set(h,'Color','r')
end
semilogx(N,-log10(err),'k-','Linewidth',2);
%print -deps2c digitsLog20000.eps