tmp_1[rb[1-cb] = u[mu][1-cb] * shift(u[nu][cb], (1-cb), FORWARD, mu);
tmp_1[rb[1-cb]] = u[mu] * shift(u[nu], FORWARD, mu);
void mciter(multi1d<LatticeColorMatrix>& u, 	    int n_over, int nheat,	    int& ntrials, int& nfails){  START_CODE("mciter");  LatticeColorMatrix tmp_1;  LatticeColorMatrix u_staple;#if 0  Real xi02;  if( AnisoP == YES )    xi02 = xi_0 * xi_0;  #else  bool AnisoP = false;  Real xi02 = 1.0;#endif  for(int iter = 0; iter <= n_over; ++iter)  {    for(int cb = 0; cb < 2; ++cb)      for(int mu = 0; mu < Nd; ++mu)      {	u_staple = 0;	for(int nu = 0; nu < Nd; ++nu)	{	  if (nu == mu) continue;	  /* Forward staple */	  /* tmp_1(x) = u(x+mu,nu)*u_dag(x+nu,mu) */	  tmp_1[rb[cb]] = shift(u[nu], FORWARD, mu) * shift(adj(u[mu]), FORWARD, nu);          if( AnisoP )  		    if( mu == t_dir || nu == t_dir )	      tmp_1[rb[cb]] *= xi02;	  /* u_staple(x) +=  tmp_1 * u_dag(x,nu)	     += u(x+mu,nu)*u_dag(x+nu,mu)*u_dag(x,nu) */	  u_staple[rb[cb]] += tmp_1 * adj(u[nu]);          	  /* Backward staple */	  /* tmp_1(x) = u(x,mu)*u(x+mu,nu) */	  tmp_1[rb[1-cb] = u[mu][1-cb] * shift(u[nu][cb], (1-cb), FORWARD, mu);          if( AnisoP == YES )  		    if( mu == t_dir || nu == t_dir )	      tmp_1[rb[1-cb]] = xi02 * tmp_1;	  /* u_staple(x) += tmp_1_dag(x-nu) * u(x-nu,nu)	     += u_dag(x+mu-nu,nu)*u_dag(x-nu,mu)*u(x-nu,nu) */          u_staple[rb[cb]] += shift(adj(tmp_1), BACKWARD, nu) * shift(u[nu], BACKWARD, nu);	}  /* closes nu loop */	if ( iter < n_over )	{	  /* Do an overrelaxation step */          /*# Loop over SU(2) subgroup index */          for(int su2_index = 0; su2_index < Nc*(Nc-1)/2; ++su2_index)            su3over(u[mu], u_staple, su2_index, rb[cb]);	}	else	{	  /* Do a heatbath step */	  /*# Loop over SU(2) subgroup index */          for(int su2_index = 0; su2_index < Nc*(Nc-1)/2; ++su2_index)	  {	    int ntry;	    int nfail;	    su3hb(u[mu], u_staple, su2_index, nheat, ntry, nfail, rb[cb]);	    ntrials = ntrials + ntry;	    nfails = nfails + nfail;	  }          	  /* Reunitarize *//*	  reunit (u[mu][cb], lbad, OPTION[REUNITARIZE_ERROR], numbad);  */	  reunit(u[mu]);	}#if 0	/* If using Schroedinger functional, reset the boundaries */	if ( SchrFun > 0 )	{	  copymask(u[mu][cb], lSFmask[mu][cb], SFBndFld[mu][cb], REPLACE);	}#endif      }        /* closes mu loop */  }    END_CODE("mciter");}