Double tgfold; Double tgfnew; Double tgf_t; Double tgf_s; Double norm; int num_sdir; bool tdirp;
LatticeColorMatrix g;
void coulGauge(multi1d<LatticeColorMatrix>& ug, int& n_gf, 	       int j_decay, const Real& GFAccu, int GFMax, 	       bool OrDo, const Real& OrPara){  Double tgfold;  Double tgfnew;  Double tgf_t;  Double tgf_s;  Double norm;  int num_sdir;  bool tdirp;  START_CODE("coulGauge");  Real xi_sq = pow(xi_0(),2);  if( j_decay >= 0 && j_decay < Nd )  {    if( tDir() >= 0 && tDir() != j_decay )    {      num_sdir = Nd - 2;      tdirp = true;      norm = Double(Layout::vol()*Nc) * (Double(num_sdir)+Double(xi_sq));    }    else    {      num_sdir = Nd - 1;      tdirp = false;      norm = Double(Layout::vol()*Nc*num_sdir);    }  }  else  {    if( tDir() >= 0 && tDir() < Nd )    {      num_sdir = Nd - 1;      tdirp = true;      norm = Double(Layout::vol()*Nc) * (Double(num_sdir)+Double(xi_sq));    }    else    {      num_sdir = Nd;      tdirp = false;      norm = Double(Layout::vol()*Nc*num_sdir);    }  }        /* Compute initial gauge fixing term: sum(trace(U_spacelike)); */  tgf_t = 0;  tgf_s = 0;  for(int mu=0; mu<Nd; ++mu)    if( mu != j_decay )    {      if( mu != tDir() )      {	tgf_s += sum(real(trace(ug[mu])));      }      else	tgf_t += sum(real(trace(ug[mu])));    }  if( tdirp )  {    tgfold = (xi_sq*tgf_t+tgf_s)/norm;    tgf_s = tgf_s/(Double(Layout::vol()*Nc*num_sdir));    tgf_t = tgf_t/(Double(Layout::vol()*Nc));  }  else  {    tgf_s = tgf_s/(Double(Layout::vol()*Nc*num_sdir));    tgfold = tgf_s;  }    n_gf = 0;  bool wrswitch = true;    /* switch for writing of gauge fixing term */  Double conver = 1;        /* convergence criterion */  multi1d<LatticeColorMatrix> u_tmp(Nd);  /* Gauge fix until converged or too many iterations */  while( toBool(conver > GFAccu)  &&  n_gf < GFMax )  {    n_gf = n_gf + 1;    if( GFMax - n_gf < 11 )       wrswitch = true;        /* Loop over checkerboards for gauge fixing */    for(int cb=0; cb<2; ++cb)    {      if (Nc > 1)      {	/* Loop over SU(2) subgroup index */	for(int su2_index=0; su2_index < Nc*(Nc-1)/2; ++su2_index)	{	  /* Gather the Nd negative links attached to a site: */	  /* u_tmp(x,mu) = U(x-mu,mu) */	  for(int mu=0; mu<Nd; ++mu)	    u_tmp[mu][rb[cb]] = shift(ug[mu], BACKWARD, mu);	  /* Now do a gauge fixing relaxation step */	  grelax(ug, u_tmp, j_decay, su2_index, cb, OrDo, OrPara);	}   /* end su2_index loop */      }      else      {	/* Gather the Nd negative links attached to a site: */	/* u_tmp(x,mu) = U(x-mu,mu) */	for(int mu=0; mu<Nd; ++mu)	  u_tmp[mu][rb[cb]] = shift(ug[mu], BACKWARD, mu);	int su2_index = -1;	/* Now do a gauge fixing relaxation step */	grelax(ug, u_tmp, j_decay, su2_index, cb, OrDo, OrPara);      }    }     /* end cb loop */    /* Reunitarize */    for(int mu=0; mu<Nd; ++mu)      reunit(ug[mu]);    /* Compute new gauge fixing term: sum(trace(U_spacelike)): */    tgf_t = 0;    tgf_s = 0;    for(int mu=0; mu<Nd; ++mu)      if( mu != j_decay )      {	if( mu != tDir() )	{	  tgf_s += sum(real(trace(ug[mu])));	}	else	  tgf_t += sum(real(trace(ug[mu])));      }    if( tdirp )    {      tgfnew = (xi_sq*tgf_t+tgf_s)/norm;      tgf_s = tgf_s/(Double(Layout::vol()*Nc*num_sdir));      tgf_t = tgf_t/(Double(Layout::vol()*Nc));    }    else    {      tgf_s = tgf_s/(Double(Layout::vol()*Nc*num_sdir));      tgfnew = tgf_s;    }    if( wrswitch )       QDPIO::cout << "COULGAUGE: iter= " << n_gf 		  << "  tgfold= " << tgfold 		  << "  tgfnew= " << tgfnew		  << "  tgf_s= " << tgf_s 		  << "  tgf_t= " << tgf_t << endl;    /* Normalized convergence criterion: */    conver = fabs((tgfnew - tgfold) / tgfnew);    tgfold = tgfnew;  }       /* end while loop */        if( wrswitch )    QDPIO::cout << "COULGAUGE: end: iter= " << n_gf 		<< "  tgfold= " << tgfold 		<< "  tgf_s= " << tgf_s 		<< "  tgf_t= " << tgf_t << endl;#if 0  /*+ debugging */  XMLBufferWriter xml_out;  push(xml_out,"Final_trace_max_in_CoulGauge");  write(xml_out, "j_decay", j_decay);  write(xml_out, "t_dir", tDir());  write(xml_out, "n_gf",n_gf);  write(xml_out, "tgfold", tgfold);  write(xml_out, "tgf_s", tgf_s);  write(xml_out, "tgf_t", tgf_t);  pop(xml_out);#endif  END_CODE("coulGauge");}