const OrderedSubset& s = M.subset();  LatticeFermion mp; Real a; Real b; Double c; Double d;  Real rsd_sq = (RsdCG * RsdCG) * Real(norm2(chi,s));       LatticeFermion r; r[s] = chi - M(M(psi, PLUS), MINUS);   LatticeFermion p; p[s] = r;   Double cp = norm2(r, s);   	       	     QDPIO::cout << "WlInvCG: k = 0  cp = " << cp << "  rsd_sq = " << rsd_sq << endl;   if ( toBool(cp  <=  rsd_sq) ) { n_count = 0; return; }     for(int k = 1; k <= MaxCG; ++k) {  c = cp;      mp[s] = M(p, PLUS);   d = norm2(mp, s);	  a = Real(c)/Real(d);   psi[s] += a * p;	     r[s] -= a * M(mp, MINUS);     cp = norm2(r, s);	                  QDPIO::cout << "WlInvCG: k = " << k << "  cp = " << cp << endl;  if ( toBool(cp  <=  rsd_sq) ) { n_count = k; return; }   b = Real(cp) / Real(c);   p[s] = r + b*p;	 } n_count = MaxCG; QDP_error_exit("too many CG iterations: count = %d", n_count);
InvCG2_a(M, chi, psi, RsdCG, MaxCG, n_count);
void InvCG2(const LinearOperator& M,	    const LatticeFermion& chi,	    LatticeFermion& psi,	    const Real& RsdCG, 	    int MaxCG, 	    int& n_count){  const OrderedSubset& s = M.subset();  LatticeFermion mp;  Real a;  Real b;  Double c;  Double d;    Real rsd_sq = (RsdCG * RsdCG) * Real(norm2(chi,s));  //                                            +  //  r[0]  :=  Chi - A . Psi[0]    where  A = M  . M      //                      +  //  r  :=  [ Chi  -  M(u)  . M(u) . psi ]  LatticeFermion r;  r[s] = chi - M(M(psi, PLUS), MINUS);  //  p[1]  :=  r[0]  LatticeFermion p;  p[s] = r;    //  Cp = |r[0]|^2  Double cp = norm2(r, s);   	       	   /* 2 Nc Ns  flops */  QDPIO::cout << "WlInvCG: k = 0  cp = " << cp << "  rsd_sq = " << rsd_sq << endl;  //  IF |r[0]| <= RsdCG |Chi| THEN RETURN;  if ( toBool(cp  <=  rsd_sq) )  {    n_count = 0;    return;  }  //  //  FOR k FROM 1 TO MaxCG DO  //  for(int k = 1; k <= MaxCG; ++k)  {    //  c  =  | r[k-1] |**2    c = cp;    //  a[k] := | r[k-1] |**2 / < p[k], Ap[k] > ;    //      	       	       	       	       	  +    //  First compute  d  =  < p, A.p >  =  < p, M . M . p >  =  < M.p, M.p >    //  Mp = M(u) * p    mp[s] = M(p, PLUS);    //  d = | mp | ** 2    d = norm2(mp, s);	/* 2 Nc Ns  flops */    a = Real(c)/Real(d);    //  Psi[k] += a[k] p[k]    psi[s] += a * p;	/* 2 Nc Ns  flops */    //  r[k] -= a[k] A . p[k] ;    //      	       +            +    //  r  =  r  -  M(u)  . Mp  =  M  . M . p  =  A . p    r[s] -= a * M(mp, MINUS);    //  IF |r[k]| <= RsdCG |Chi| THEN RETURN;    //  cp  =  | r[k] |**2    cp = norm2(r, s);	                /* 2 Nc Ns  flops */    QDPIO::cout << "WlInvCG: k = " << k << "  cp = " << cp << endl;    if ( toBool(cp  <=  rsd_sq) )    {      n_count = k;      return;    }    //  b[k+1] := |r[k]|**2 / |r[k-1]|**2    b = Real(cp) / Real(c);    //  p[k+1] := r[k] + b[k+1] p[k]    p[s] = r + b*p;	/* Nc Ns  flops */  }  n_count = MaxCG;  QDP_error_exit("too many CG iterations: count = %d", n_count);}