% Generates Table 7.5 

format short g

n = 2000000; p = 32452843; 

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;

table = []; kmax = 28;
for n = 20000*[1 10 100]
    A = AA(1:n,1:n); b = bb(1:n);
    D=spdiags(diag(A),0,n,n);
    for tol = [1e-11,1.5e-15]
        [x,fail,r,k] = pcg(A,b,tol,kmax,D);
        if not(fail), 
            x1 = x(1); 
            err = r/x1;
            m = 2*floor(log2(n))+4; u = eps/2; gamma = m*u/(1-m*u);
            err2 = (norm(r)+gamma*norm(A*abs(x)+abs(b)))/x(1);
            table = [table; n tol k err err2];
        end
    end
end

table