function x=cgm(A,b)

%%% Solves Ax=b using the conjugate gradient method.
%%% The matrix A should be symmetric and positive definite.

%%% Usage (examples at end of file):
%   x=cgm(A,b);

% Required arguments:
%  A: the nxn matrix
%  b: column n-vector

%%% tol = error tolerance; can be adjusted as needed
tol=2*eps;
%tol=1e-12;

%%% Imax = max number of iteration steps allowed
n=length(b);
Imax=2*n;

% the cgm code begins
nm=length(b);
x=zeros(nm,1);
r=b;
d=r;
rr=r'*r;
err=3*tol;
nb=norm(b,inf);
if nb==0
    error('b=0 to machine precision')
end
beta=0;
iter=0;
while err>tol
    d=r+beta*d;
    q=A*d;
    dq=d'*q;
    alpha=rr/dq;
    ad=alpha*d;
    x=x+ad;
    r=r-alpha*q;
    rr0=rr;
    rr=r'*r;
    err=norm(ad,inf);
    %err=norm(r,inf)/nb;
    beta=rr/rr0;
    iter=iter+1;
    if iter>Imax
        error('Method failed: exceeded maximum number of iterations')
    end
end
%fprintf('\n  Number of CGM Iterations =  %i \n\n',iter)

%%% Example 1: random SPD matrix
%
% n=1000;
% A1=rand(n,n);             %% generate A
% A=A1'+A1+diag(2*n*ones(n,1));
% xe=2*(1-0.5*rand(n,1));     %% pick solution
% b=A*xe;
% x=cgm(A,b);
% error=norm(x-xe,inf)/norm(xe,inf)  %% check accuracy
%
%%%% The error in the above is 6e-15

%%% Example 2: arrowhead matrix
%
% n=1000;                   %% generate A
% A=diag(2*ones(n,1));
% A(2:n,1)=1; A(1,2:n)=1; A(1,1)=n+1;
% xe=2*(1-0.5*rand(n,1));     %% pick solution
% b=A*xe;
% x=cgm(A,b);
% error=norm(x-xe,inf)/norm(xe,inf)  %% check accuracy
% %
%%%% The error in the above is 8e-16

%  Background
%  The algorithm is described in the text "Introduction to Scientific
%  Computing and Data Analysis" (Holmes, 2016).

%  Version history
%  verion 1.0: March 28, 2016
%  verion 1.1: April 15, 2018 (modified the examples, CGM code unchanged)