function x=cgmS(C,b)

%%% Solves Ax=b using the conjugate gradient method; it's capable of   
%%% quickly solving equations with very large sparse matrices (even those  
%%% that exceed MATLAB's storage limitations).  The matrix A should be  
%%% symmetric and positive definite.  

%%% usage (examples at end of file):
%   x=cgmS(C,b);

% Required arguments:
%  C: this is the matrix A stored in a row-compressed format 
%     (this is explained at end of file)
%  b: column n-vector

%%% Comment: This program is most useful for large sparse matrices.  For 
%%% non-large matrices it might be easier, and faster, to use cgm(A,b), which
%%% is also available via MATLAB's file exchange.

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

%%% Imax = max number of iteration steps allowed
n=length(b);
Imax=n+10;

v(1)=0; 
for i=1:n
    v(i+1)=find(C(2,:)==i,1,'last');
    c{i}=C(1,v(i)+1:v(i+1));
    ddC{i}=C(3,v(i)+1:v(i+1));
end
x=zeros(n,1);
r=b;
d=r;
rr=r'*r;
beta=0;
err=100*tol+100;
nb=norm(b,inf);
if nb==0
    return
end
q=zeros(n,1);
iter=0;
while err>tol 
    d=r+beta*d;
    for i=1:n
        q(i)=c{i}*d(ddC{i});
    end
    dq=d'*q;
    alpha=rr/dq;
    x=x+alpha*d;
    r=r-alpha*q;
    rr0=rr;
    rr=r'*r;
    %err=norm(alpha*d,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

%%% Example 1: random sparse matrix
%
% n=10000;
% A1=zeros(n,n);             %% generate A
% for i=1:n
%     A1(i,i)=2+randi(8);
%     A1(randi(n),randi(n))=1;
% end
% A=A1'*A1;
% [row,col]=find(A~=0);      %% find C from A
% CR=[row col];
% RC=sortrows(CR);
% for i=1:length(row)
%     C(1,i)=A(RC(i,1),RC(i,2));
%     C(2,i)=RC(i,1);
%     C(3,i)=RC(i,2);
% end
% xe=2*(1-0.5*rand(n,1));     %% pick solution
% b=A*xe;
% x=cgmS(C,b);
% error=norm(x-xe,inf)/norm(xe,inf)  %% check accuracy
%
%%%% The error in the above is 7.9e-15

%%% Example 2: arrowhead matrix
%
% n=10000;                   %% generate A
% A=zeros(n,n);
% for i=2:n-1
%     A(i,i)=2;
%     A(1,i)=1;
%     A(i,1)=1;
% end
% A(1,1)=2*n;
% A(n,n)=2*n;
% [row,col]=find(A~=0);      %% find C from A
% CR=[row col];
% RC=sortrows(CR);
% for i=1:length(row)
%     C(1,i)=A(RC(i,1),RC(i,2));
%     C(2,i)=RC(i,1);
%     C(3,i)=RC(i,2);
% end
% b=2*(1-0.5*rand(n,1));     %% pick b
% x=cgmS(C,b);
% xx=A\b;                    %% compare with MATLAB solution
% error=norm(x-xx,inf)/norm(xx,inf)
%
%%%% The error in the above is 8e-15

%%% Format of Matrix C
%%% This is 3xN, where N=total number of nonzero entries in A
%%% first row: all of the nonzero entries in original matrix A 
%%% 2nd and 3d rows: contain the i (2nd row) and j (3d row) coordinates from  
%%%     A for the nonzero entry in the first row of C
%%% Example:
%%% A=[1 0 0; 0 3 2; 0 2 5];
%%% C(1,:)=[1 3 2 2 5];
%%% C(2,:)=[1 2 2 3 3];
%%% C(3,:)=[1 2 3 2 3];
%%% A procedure for computing C from A is given in the above two examples 
%%% (lines 75-82 and 101-108).  However for very large matrices (larger than 
%%% MATLAB normally works with), you will need to produce C directly.

%  Background
%  This was written to solve large sparse matrix equations. The 
%  algorithm is described in the text "Introduction to Scientific Computing 
%  and Data Analysis" (Holmes, 2016).  

%  It has been tested using MATLAB, version R2016a
%  March 28, 2016
















