The conjugate gradient algorithm

To solve the linear system AU=F we use the conjugate gradient algorithm:

initialization: U=0, H = G  = F, eps = 1e-6  
iterations:     d = |G|
                r = G.AH / H.AH
                U := U + r H
                G := G + r
                c  = |G|/d
                H := G + c H
                if(|G|<eps |F|) stop

Here we notice that

where , so that this is programmed as

Scalar products and norms are programmed as usual

float scal = 0, norm = 0;
for(int i=0; i<g.nv; i++)  
{  scal += a[i] * b[i];  // scalar product a * b
   norm += pow(a[i],2);  // norm of a
}
norm = sqrt(norm);