*/
void MInvCG_a(const LinearOperator<T>& A, 	      const T& chi, 	      multi1d<T>& psi,	      const multi1d<Real>& shifts, 	      const multi1d<Real>& RsdCG, 	      int MaxCG,	      int& n_count){  START_CODE("MinvCG");  const OrderedSubset& sub = A.subset();  int n_shift = shifts.size();  if (n_shift == 0) {    QDP_error_exit("MinvCG: You must supply at least 1 mass: mass.size() = %d",		   n_shift);  }  /* Now find the smallest mass */  int isz = 0;  for(int findit=1; findit < n_shift; ++findit) {    if ( toBool( shifts[findit] < shifts[isz])  ) {       isz = findit;    }  }#if 0   QDPIO::cout << "n_shift = " << n_shift << " isz = " << isz << " shift = " << shifts[0] << endl;#endif  // We need to make sure, that psi is at least as big as the number  // of shifts. We resize it if it is not big enough.  // However, it is allowed to be bigger.  if( psi.size() <  n_shift ) {       psi.resize(n_shift);  }  // For this algorithm, all the psi have to be 0 to start  // Only is that way the initial residuum r = chi  for(int i= 0; i < n_shift; ++i) {     psi[i][sub] = zero;  }    // If chi has zero norm then the result is zero  Double chi_norm = sqrt(norm2(chi,sub));  if( toBool( chi_norm < fuzz )) {     n_count = 0;    // The psi are all zero anyway at this point    // for(int i=0; i < n_shift; i++) { psi[i] = zero; }    return;  }  multi1d<Double> rsd_sq(n_shift);  multi1d<Double> rsdcg_sq(n_shift);  Double cp = chi_norm * chi_norm;  int s;  for(s = 0; s < n_shift; ++s)  {    rsdcg_sq[s] = RsdCG[s] * RsdCG[s];  // RsdCG^2    rsd_sq[s] = Real(cp) * rsdcg_sq[s]; // || chi ||^2 RsdCG^2  }    // r[0] := p[0] := Chi   T r;  r[sub] = chi;  // Psi[0] := 0;  multi1d<T> p(n_shift);  for(s = 0; s < n_shift; ++s) {    p[s][sub] = chi;  }  //  b[0] := - | r[0] |**2 / < p[0], Ap[0] > ;/  //  First compute  d  =  < p, A.p >   //  Ap = A . p  */  LatticeFermion Ap;  A(Ap, p[isz], PLUS);  Ap[sub] += p[isz] * shifts[isz];  /*  d =  < p, A.p >  */  Double d = real(innerProduct(p[isz], Ap, sub)); // 2Nc Ns flops      Double b = -cp/d;  /* Compute the shifted bs and z */  multi1d<Double> bs(n_shift);  multi2d<Double> z(2, n_shift);  int iz;  z[0][isz] = Double(1);  z[1][isz] = Double(1);  bs[isz] = b;  iz = 1;  for(s = 0; s < n_shift; ++s)  {    if( s != isz ) {      z[1-iz][s] = Double(1);      z[iz][s] = Double(1) / (Double(1) - (Double(shifts[s])-Double(shifts[isz]))*b);      bs[s] = b * z[iz][s];    }  }  //  r[1] += b[0] A . p[0];   r[sub] += Ap * Real(b);	                        // 2 Nc Ns  flops  //  Psi[1] -= b[0] p[0] = - b[0] chi;  for(s = 0; s < n_shift; ++s) {    psi[s][sub] = - Real(bs[s])*chi;                      //  2 Nc Ns  flops   }    //  c = |r[1]|^2     Double c = norm2(r,sub);   	       	         //  2 Nc Ns  flops   // Check convergence of first solution  multi1d<bool> convsP(n_shift);  for(s = 0; s < n_shift; ++s) {    convsP[s] = false;  }  bool convP = toBool( c < rsd_sq[isz] );#if 0   QDPIO::cout << "MInvCG: k = 0  r = " << sqrt(c) << endl;#endif  //  FOR k FROM 1 TO MaxCG DO  //  IF |psi[k+1] - psi[k]| <= RsdCG |psi[k+1]| THEN RETURN;   Double z0, z1;  Double ztmp;  Double cs;  Double a;  Double as;  Double  bp;  int k;    for(k = 1; k <= MaxCG && !convP ; ++k)  {    //  a[k+1] := |r[k]|**2 / |r[k-1]|**2 ;     a = c/cp;    //  p[k+1] := r[k+1] + a[k+1] p[k];     //  Compute the shifted as */    //  ps[k+1] := zs[k+1] r[k+1] + a[k+1] ps[k];    for(s = 0; s < n_shift; ++s) {      if (!convsP[s]) {		if (s == isz) {	  p[s][sub] *= Real(a);	                              // Nc Ns  flops 	  p[s][sub] += r;	                              // Nc Ns  flops 	}	else {	  as = a * z[iz][s]*bs[s] / (z[1-iz][s]*b);	  	  p[s][sub] *= Real(as);	                             // Nc Ns  flops 	  p[s][sub] += r * Real(z[iz][s]);	                     // Nc Ns  flops 	}      }    }    //  cp  =  | r[k] |**2     cp = c;    //  b[k] := | r[k] |**2 / < p[k], Ap[k] > ;    //  First compute  d  =  < p, A.p >      //  Ap = A . p     A(Ap, p[isz], PLUS);    Ap[sub] += p[isz] * shifts[isz];    /*  d =  < p, A.p >  */    d = real(innerProduct(p[isz], Ap, sub));                   //  2 Nc Ns  flops        bp = b;    b = -cp/d;    // Compute the shifted bs and z     bs[isz] = b;    iz = 1 - iz;    for(s = 0; s < n_shift; s++) {                  if (s != isz && !convsP[s] ) {	z0 = z[1-iz][s];	z1 = z[iz][s];	z[iz][s] = z0*z1*bp;	z[iz][s] /= b*a*(z1-z0) + z1*bp*(Double(1) - (shifts[s] - shifts[isz])*b);	bs[s] = b*z[iz][s]/z0;      }    }    //  r[k+1] += b[k] A . p[k] ;     r[sub] += Ap * Real(b);	        // 2 Nc Ns  flops    //  Psi[k+1] -= b[k] p[k] ;     for(s = 0; s < n_shift; ++s) {      if (! convsP[s] ) {	psi[s][sub] -= p[s] * Real(bs[s]);	// 2 Nc Ns  flops       }    }    //  c  =  | r[k] |**2     c = norm2(r,sub);	                // 2 Nc Ns  flops     //    IF |psi[k+1] - psi[k]| <= RsdCG |psi[k+1]| THEN RETURN;    // or IF |r[k+1]| <= RsdCG |chi| THEN RETURN;    convP = true;    for(s = 0; s < n_shift; s++) {      if (! convsP[s] ) {	// Convergence methods #if 0	// Check norm of shifted residuals 	Double css = c * z[iz][s]* z[iz][s];#if 0		if (s == isz) {	  QDPIO::cout << "MInvCG: k = " << k <<"  r =  " 		      << sqrt(css/chi_norm) << endl;	}#endif 	convsP[s] = toBool(  css < rsd_sq[s] );#else	// Check relative error of solution 	cs = norm2(p[s],sub);         	        // 2 Nc Ns  flops 	d = norm2(psi[s],sub);         	        // 2 Nc Ns  flops 	cs *= bs[s]*bs[s];	d *= rsdcg_sq[s];	convsP[s] = toBool( cs < d );#if 0	if (s == isz) {	  QDPIO::cout  << "MInvCG: k = " << k << " s = " << s << " r = " 		       << sqrt(cs) << " d = " << sqrt(d) << endl;	}#endif#endif      }      convP &= convsP[s];    }    n_count = k;  }#if 1  // Expicitly check the solutions  for(s = 0; s < n_shift; ++s)  {    A(Ap, psi[s], PLUS);    Ap[sub] += psi[s] * shifts[s];    Ap[sub] -= chi;    c = norm2(Ap,sub);	                /* 2 Nc Ns  flops */    QDPIO::cout << "MInvCG (conv): s = " << s                 << " shift = " << shifts[s]		<< " r = " <<  Real(sqrt(c/chi_norm)) << endl;		  }  /* end */#endif  if (n_count == MaxCG) {    QDP_error_exit("too many CG iterationns: %d\n", n_count);  }  END_CODE("MinvCG");  return;  }