function [x,k]=CGr(A,b,x0)
% Gradientes Conjugados para resolver Ax=b
%
 if nargin<3, x0=zeros(length(b),1); end
 k=0; x=x0; r=b-A*x; d=r'*r; p=r; tol=sqrt(eps)*d;
 while d>tol
    k=k+1;
    v=A*p;
    alpha=d/(p'*v);
    x=x+alpha*p;
    r=r-alpha*v;
    beta=d;
    d=r'*r;
    beta=d/beta;
    p=r+beta*p;
 end
end