%--------------- function [err_sol,resvec] = mycg(A,b,x,iter) %input: matrix SPD A, rhs b, initial vector x %output: sequence of vectors with errors (in energy norm) and norms of residuals % sol = A\b; r = b - A*x; p = r; resvec = zeros(iter,1); resvec(1) = norm(r); for j=2:iter a = (r'*r)/(p'*A*p); x = x+a*p; z= r'*r; r = r - a*A*p; b = (r'*r)/z; p = r+b*p; resvec(j) = norm(r); err_sol(j) = sqrt((sol-x)'*A*(sol-x)); end; return end