%
% Conjugate gradient method for Ax=b with stopping criterion
% based on the eigenvalue estimation (condition number) of A.
%
function x=cgeig(A,b,x0,toler,itmax)
x=x0;
bnorm = norm(b);
r = b-A*x;
p = r;
rr = r'*r;
kappa = 1e6;
T = sparse(itmax+1,itmax+1);   
for iter=1:itmax
    Ap = A*p;
    alpha = rr/(p'*Ap);
    x = x + alpha*p;
    r = r - alpha*Ap;
    rrnew = r'*r;
    if sqrt(rrnew) < (toler/kappa)*bnorm
        break;
    end
    beta = rrnew/rr;
    p = r + beta*p;
    rr = rrnew;
    %
    if iter==1
        T(iter,iter) = 1/alpha;
    else    
        T(iter,iter) = 1/alpha + bold/aold;
    end
    T(iter,iter+1) = sqrt(beta)/alpha;
    T(iter+1,iter) = T(iter,iter+1);
    
    if (iter>1 && iter <=10) || mod(iter,20)==0
        %
        % Compute smallest and largest eigenvalue of the tridiagonal matrix
        % T (this is cheap). (See: Saad page 203)
        % 
        %
        lam = eigs(T(1:iter,1:iter),2,'BE');  
        kappa = max(lam)/min(lam);
    end
    aold = alpha;
    bold = beta;
end
fprintf(1,'Iterations = %d. Condition estimate = %e\n',iter,kappa);
