function [lamb, y, res, info]=CG(A, x, k, U, nrmA, tol)
%
%  Conjugate Gradient method to compute the kth smallest eigenvalue of A
%
%
% RCL 6/28/2009, updated 7/31/2009
%---------------------------------------------------------------
%
%  Input
%
%        A     n-by-n Hermitian matrix
%        x     n-vector, initial guess to the kth eigenvector
%        k     index of the desired eigenvalue
%        U     n-by-(k-1) the 1st (k-1) accurate eigenvectors
%        nrmA  An estimate of ||A||
%        tol   tolerance for testing convergence
%
%  Output
%
%        lamb  converged eigenvalue
%        y     corresponding eigenvector
%        res   residual error  ||A y - lamb y||
%        info  the number of CG steps
%
% Note: For the standard symmetric eigenvalue problem, if x'*U=0, then
% r'*U=0, too. But numberically, since x'*U approx 0, it is possible
% after many SD steps x may have nonnegligible commonents in U. Therefore
% it may be necessary to orthogonalize x against U every many SD steps.
%
n=max(size(A));
maxitn=round(0.1*n); whichbeta = 1; 
for i=1:k-1,
    s=U(:,i)'*x; x=x-s*U(:,i);
end
x0=(1/norm(x))*x; Ax=A*x0;
rho0=x0'*Ax; r0=Ax-rho0*x0; nrmr0=norm(r0);
err=nrmr0/(nrmA+abs(rho0)); itn=0;
p0=r0; 

whichbeta=1;

while err > tol & itn < maxitn,
    for i=1:k-1,
        s=U(:,i)'*p0; p0=p0-s*U(:,i);
    end
    nrmp0=norm(p0); p0n=(1/nrmp0)*p0;
    rhos=[rho0 p0n'*(A*p0n)];
    [x1, t]=onestep0(Ax,x0,p0n,rhos);
    for i=1:k-1,
        s=U(:,i)'*x1; x1=x1-s*U(:,i);
    end
    nrmx1=norm(x1); x1=(1/nrmx1)*x1; 
    Ax=A*x1; rho1=x1'*Ax; r1=Ax-rho1*x1;
    if whichbeta == 1,
       beta=r1'*r1/(nrmr0*nrmr0); 
    else
       tmp=nrmx1*r1;
       beta=r1'*(tmp-r0)/(nrmr0*nrmr0); 
    end
    p1=r1+beta*p0;
    x0=x1; rho0=rho1; r0=r1; p0=p1; 
    nrmr0=norm(r0); err=nrmr0/(nrmA+abs(rho0));
    itn = itn+1;
end
info = itn; y=x0; lamb=rho0; res=nrmr0;


