const Real& rho,
void InvRelCG2_a(const ApproxLinearOperator<T>& M,		 const T& chi,		 T& psi,		 const Real& rho,		 const Real& RsdCG, 		 int MaxCG, 		 int& n_count){  const OrderedSubset& s = M.subset();    //  Real rsd_sq = (RsdCG * RsdCG) * Real(norm2(chi,s));  Real chi_sq =  Real(norm2(chi,s));  QDPIO::cout << "chi_norm = " << sqrt(chi_sq) << endl;  Real rsd_sq = (RsdCG * RsdCG) * chi_sq;  //                                            +  //  r[0]  :=  Chi - A . Psi[0]    where  A = M  . M      //                      +  //  r  :=  [ Chi  -  M(u)  . M(u) . psi ]  T  r, mp, mmp;  r=zero;  mp=zero;  mmp=zero;  M(mp, psi, PLUS);  M(mmp, mp, MINUS);  r[s] = chi - mmp;  //  p[1]  :=  r[0]  T p;  p = zero;  p[s] = r;    //  Cp = |r[0]|^2  Double c = norm2(r, s);   	       	   /* 2 Nc Ns  flops */  Double cp = c;  Double zeta = Double(1)/c;  QDPIO::cout << "InvRelCG2: k = 0  c = " << cp << "  rsd_sq = " << rsd_sq << endl;  //  IF |r[0]| <= RsdCG |Chi| THEN RETURN;  if ( toBool(c  <=  rsd_sq) )  {    n_count = 0;    return;  }  //  //  FOR k FROM 1 TO MaxCG DO  //  Real a;  Double d;    for(int k = 1; k <= MaxCG; ++k)  {    // Inner tolerance = epsilon || chi || || p || sqrt(zeta) / 2    Real inner_tol = RsdCG*sqrt(chi_sq)*sqrt(norm2(p))*sqrt(zeta)/(Real(8)*rho);    // Compute M^{dag} M p    M(mp, p, PLUS, inner_tol);    M(mmp, mp, MINUS, inner_tol);    //  d = < M^{dag} M p, p>    d = innerProductReal(mmp, p, s);	/* 2 Nc Ns  flops */    //  a[k] := | r |**2 / < p[k], Ap[k] > ;    a = Real(c)/Real(d);    //  Psi[k] += a[k] p[k]    psi[s] += a * p;	/* 2 Nc Ns  flops */    //  r[k] -= a[k] M^{dag}M. p[k] ;    r[s] -= a * mmp;    //  IF |r[k]| <= RsdCG |Chi| THEN RETURN;    //  cp  =  | r[k] |**2    c = norm2(r, s);	                /* 2 Nc Ns  flops */    // update relaxation factor    zeta += Double(1)/c;    // Update p as:    // p[k+1] := r[k] + c/cp * p    //    // we put c/cp into b    //    //  b[k+1] := |r[k]|**2 / |r[k-1]|**2    Real b = Real(c) / Real(cp);    //  p[k+1] := r[k] + b[k+1] p[k]    p[s] = r + b*p;	/* Nc Ns  flops */    cp = c;    QDPIO::cout << "InvCG: k = " << k << "  cp = " << cp << endl;    if ( toBool(cp  <=  rsd_sq) )    {      n_count = k;      return;    }  }  n_count = MaxCG;  QDP_error_exit("too many CG iterations: count = %d", n_count);}