push(xml_out, "time_direction"); write(xml_out, "t_dir",t_dir); pop(xml_out);
void dwf_quarkProp4_a(LatticePropagator& q_sol, 		      XMLWriter& xml_out,		      const LatticePropagator& q_src,		      const int t_src,		      const C<T>& S_f,		      Handle<const ConnectState> state,		      enum InvType invType,		      const Real& RsdCG, 		      int MaxCG, int& ncg_had){  START_CODE("dwf_quarkProp4");  push(xml_out, "DWF_QuarkProp4");  ncg_had = 0;  multi1d<LatticePropagator> prop5d(S_f.size()) ;  LatticePropagator q_mp  ;  multi1d<LatticeFermion> psi(S_f.size()) ;      // This version loops over all color and spin indices  for(int color_source = 0; color_source < Nc; ++color_source)  {    for(int spin_source = 0; spin_source < Ns; ++spin_source)    {      QDPIO::cout<<"dwf_quarkProp4:: doing color  : "<< color_source;      QDPIO::cout<<" and spin : "<< spin_source<<endl  ;      psi = zero ;  // note this is ``zero'' and not 0      LatticeFermion tmp ;      tmp = zero ;      PropToFerm(q_src, tmp, color_source, spin_source);	         /*        * Normalize the source in case it is really huge or small -        * a trick to avoid overflows or underflows       */      Real fact = Real(1) / sqrt(norm2(tmp));      tmp *= fact ;      //QDPIO::cout<<"Normalization Factor: "<< fact<<endl ;      int N5(S_f.size());      multi1d<LatticeFermion> chi(N5) ;      chi = zero ;      // Split the source to oposite walls according to chirality      chi[0   ] = chiralProjectPlus(tmp) ;      chi[N5-1] = chiralProjectMinus(tmp) ;       // now we are ready invert	         int n_count;      // Compute the propagator for given source color/spin.	         S_f.qpropT(psi, state, chi, invType, RsdCG, MaxCG, n_count);      ncg_had += n_count;      push(xml_out,"Qprop");      write(xml_out, "RsdCG", RsdCG);      Write(xml_out, n_count);      pop(xml_out);      // Unnormalize the source following the inverse       // of the normalization above      fact = Real(1) / fact;      for(int i = 0; i < N5; ++i)	psi[i] *= fact;       /*       * Move the solution to the appropriate components       * of quark propagator.       */	         //First the 5D quark propagator      for(int s(0);s<S_f.size();s++)	FermToProp(psi[s], prop5d[s], color_source, spin_source);      // Now get the 4D propagator too      tmp = chiralProjectMinus(psi[0]) + chiralProjectPlus(psi[N5-1]) ;      // move solution to the appropriate components of the 4d      // quark propagator      FermToProp(tmp, q_sol, color_source, spin_source);      // move solution to the appropriate components of the 4d      // midpoint quark propagator       tmp = chiralProjectPlus(psi[N5/2 - 1]) + chiralProjectMinus(psi[N5/2]) ;      FermToProp(tmp, q_mp, color_source, spin_source);	       }	/* end loop over spin_source */  } /* end loop over color_source */  // constuct the conserved axial current correlator  LatticeComplex cfield ;  dwf_conserved_axial_ps_corr(cfield,state->getLinks(),prop5d,3);			         multi1d<DComplex> corr ;       SftMom trick(0,false,3) ;     corr = sumMulti(cfield, trick.getSet());  // Length of lattice in decay direction  int length = trick.numSubsets();  multi1d<Real> mesprop(length);  //multi1d<Complex> axps(length);  for(int t(0);t<length; t++){    int t_eff( (t - t_src + length) % length ) ;    mesprop[t_eff] = real(corr[t]) ;   }  //push(xml_out, "t_dir");  //Write(xml_out, 3);  //pop(xml_out);  push(xml_out, "DWF_ConservedAxial");  //write(xml_out, "corr", corr);  Write(xml_out, mesprop);   pop(xml_out);  // The local axial corruent pseudoscalar correlator  cfield = trace( adj(q_sol)*Gamma(8)*q_sol ) ;  corr = sumMulti(cfield, trick.getSet()) ;  for(int t(0);t<length; t++){    int t_eff( (t - t_src + length) % length ) ;    mesprop[t_eff] = real(corr[t]) ;   }  push(xml_out, "DWF_LocalAxial");  Write(xml_out, mesprop);   pop(xml_out);  //Now the midpoint Pseudoscalar correlator  multi1d<Double> tmp(length);  tmp = sumMulti(localNorm2(q_mp), trick.getSet());  for(int t(0);t<length; t++){    int t_eff( (t - t_src + length) % length ) ;    mesprop[t_eff] = tmp[t] ;   }    push(xml_out, "DWF_MidPoint_Pseudo");  Write(xml_out, mesprop);  pop(xml_out);  tmp = sumMulti(localNorm2(q_sol), trick.getSet());  for(int t(0);t<length; t++){    int t_eff( (t - t_src + length) % length ) ;    mesprop[t_eff] = tmp[t] ; // only need the zero momentum  }  push(xml_out, "DWF_Psuedo_Pseudo");  Write(xml_out, mesprop);  pop(xml_out);  pop(xml_out);   // DWF_QuarkProp4  check_dwf_ward_identity(state->getLinks(),prop5d,q_sol,q_mp,S_f.quark_mass());  END_CODE("dwf_quarkProp4");}