% The standard Conjugate Gradient implementation
% with a matrix A. No preconditioning.
% Parameters: zero initial guess; b-right_hand_side; y-the solution
% ---------------------------------------------------------------
function [it,y,resvec]=cg(A,b,eps_cg)
it=1;
y=zeros(size(A,1),1);
resvec=[];
g=A*y-b;
qq0=g'*g;
disp(['___________cg_delta0:',num2str(qq0)])
ceps = eps_cg*qq0;
qq1=1e+12;
resvec(it)=qq0;
d=-g;
    while qq1>ceps,
        h=A*d;
        tau=qq0/(d'*h);
        y=y+tau*d,tau
        g=g+tau*h;
        qq1=g'*g;
        disp(['______________________________cg_delta1:',num2str(qq1)])
        beta=qq1/qq0;
        d=-g+beta*d;
        it=it+1;
        qq0=qq1;
	resvec(it)=qq0;
    end
it=it-1;
return
% ---------------------------------------------------------------