diff --git a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h index 2300afd3..27931f91 100644 --- a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h +++ b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h @@ -60,6 +60,50 @@ public: // virtual void Instantiatable(void)=0; virtual void Instantiatable(void) =0; + void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector boundary, std::vector twist) + { + std::cout << "Free Propagator for PartialFraction"<_fdimensions[nu+shift]))); + //momenta for propagator shifted by twist+boundary + twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI)); + } + in_buf = exp(ci*ph*(-1.0))*in; + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + + //phase for boundary condition + out = out * exp(ci*ph); + }; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + std::vector twist(Nd,0.0); //default: periodic boundarys in all directions + std::vector boundary; + for(int i=0;i &out); diff --git a/Grid/qcd/action/fermion/PartialFractionFermion5D.h b/Grid/qcd/action/fermion/PartialFractionFermion5D.h index 54f8547f..34aa80b8 100644 --- a/Grid/qcd/action/fermion/PartialFractionFermion5D.h +++ b/Grid/qcd/action/fermion/PartialFractionFermion5D.h @@ -39,7 +39,7 @@ class PartialFractionFermion5D : public WilsonFermion5D public: INHERIT_IMPL_TYPES(Impl); - const int part_frac_chroma_convention=1; + const int part_frac_chroma_convention=0; void Meooe_internal(const FermionField &in, FermionField &out,int dag); void Mooee_internal(const FermionField &in, FermionField &out,int dag); @@ -83,12 +83,63 @@ public: GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD M5,const ImplParams &p= ImplParams()); + PartialFractionFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5,std::vector &_qmu,const ImplParams &p= ImplParams()); + + void FreePropagator(const FermionField &in,FermionField &out,RealD mass,std::vector boundary, std::vector twist) + { + std::cout << "Free Propagator for PartialFraction"<_fdimensions[nu+shift]))); + //momenta for propagator shifted by twist+boundary + twist[nu] = twist[nu] + boundary_phase/((2.0*M_PI)); + } + in_buf = exp(ci*ph*(-1.0))*in; + + theFFT.FFT_all_dim(in_k,in,FFT::forward); + this->MomentumSpacePropagatorHw(prop_k,in_k,mass,twist); + theFFT.FFT_all_dim(out,prop_k,FFT::backward); + + //phase for boundary condition + out = out * exp(ci*ph); + }; + + virtual void FreePropagator(const FermionField &in,FermionField &out,RealD mass) { + std::vector twist(Nd,0.0); //default: periodic boundarys in all directions + std::vector boundary; + for(int i=0;i qmu; RealD mass; RealD dw_diag; RealD R; diff --git a/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h index 6687800e..4bfbd31e 100644 --- a/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/ContinuedFractionFermion5DImplementation.h @@ -42,13 +42,13 @@ template void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata) { // How to check Ls matches?? - // std::cout<n << " - n"<da << " -da "<db << " -db"<dn << " -dn"<dd << " -dd"<n << " - n"<da << " -da "<db << " -db"<dn << " -dn"<dd << " -dd"<Ls; + std::cout<db==Ls);// Beta has Ls coeffs R=(1+this->mass)/(1-this->mass); @@ -320,7 +320,7 @@ ContinuedFractionFermion5D::ContinuedFractionFermion5D( int Ls = this->Ls; conformable(solution5d.Grid(),this->FermionGrid()); conformable(exported4d.Grid(),this->GaugeGrid()); - ExtractSlice(exported4d, solution5d, Ls-1, Ls-1); + ExtractSlice(exported4d, solution5d, Ls-1, 0); } template void ContinuedFractionFermion5D::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d) @@ -330,7 +330,7 @@ ContinuedFractionFermion5D::ContinuedFractionFermion5D( conformable(input4d.Grid() ,this->GaugeGrid()); FermionField tmp(this->FermionGrid()); tmp=Zero(); - InsertSlice(input4d, tmp, Ls-1, Ls-1); + InsertSlice(input4d, tmp, Ls-1, 0); tmp=Gamma(Gamma::Algebra::Gamma5)*tmp; this->Dminus(tmp,imported5d); } diff --git a/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h index 0206828b..cbbda785 100644 --- a/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h @@ -255,15 +255,76 @@ void PartialFractionFermion5D::M_internal(const FermionField &psi, Fermi } { + // The 'conventional' Cayley overlap operator is + // + // Dov = (1+m)/2 + (1-m)/2 g5 sgn Hw + // + // + // With massless limit 1/2(1+g5 sgnHw) + // + // Luscher shows quite neatly that 1+g5 sgn Hw has tree level propagator i qslash +O(a^2) + // + // However, the conventional normalisation has both a leading order factor of 2 in Zq + // at tree level AND a mass dependent (1-m) that are convenient to absorb. + // + // In WilsonFermion5DImplementation.h, the tree level propagator for Hw is + // + // num = -i sin kmu gmu + // + // denom ( sqrt(sk^2 + (2shk^2 - 1)^2 + // b_k = sk2 - M5; + // + // w_k = sqrt(sk + b_k*b_k); + // + // denom= ( w_k + b_k + mass*mass) ; + // + // denom= one/denom; + // out = num*denom; + // + // Chroma, and Grid define partial fraction via 4d operator + // + // Dpf = 2/(1-m) x Dov = (1+m)/(1-m) + g5 sgn Hw + // + // Now since: + // + // (1+m)/(1-m) = (1-m)/(1-m) + 2m/(1-m) = 1 + 2m/(1-m) + // + // This corresponds to a modified mass parameter + // + // It has an annoying + // + // double R=(1+this->mass)/(1-this->mass); //R g5 psi[Ls] + p[0] H ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1); - + for(int b=0;b::SetCoefficientsZolotarev(RealD zolo_hi,App int Ls = this->Ls; conformable(solution5d.Grid(),this->FermionGrid()); conformable(exported4d.Grid(),this->GaugeGrid()); - ExtractSlice(exported4d, solution5d, Ls-1, Ls-1); + ExtractSlice(exported4d, solution5d, Ls-1, 0); } template void PartialFractionFermion5D::ImportPhysicalFermionSource(const FermionField &input4d,FermionField &imported5d) @@ -421,7 +482,8 @@ void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,App conformable(input4d.Grid() ,this->GaugeGrid()); FermionField tmp(this->FermionGrid()); tmp=Zero(); - InsertSlice(input4d, tmp, Ls-1, Ls-1); + std::cout << " importing to slice " << Ls-1 <Dminus(tmp,imported5d); } @@ -442,7 +504,7 @@ PartialFractionFermion5D::PartialFractionFermion5D(GaugeField &_Umu, { int Ls = this->Ls; - + qmu.resize(0); assert((Ls&0x1)==1); // Odd Ls required int nrational=Ls-1; @@ -460,6 +522,22 @@ PartialFractionFermion5D::PartialFractionFermion5D(GaugeField &_Umu, Approx::zolotarev_free(zdata); } +template +PartialFractionFermion5D::PartialFractionFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5, + std::vector &_qmu, + const ImplParams &p) + : PartialFractionFermion5D(_Umu, + FiveDimGrid,FiveDimRedBlackGrid, + FourDimGrid,FourDimRedBlackGrid, + _mass,M5,p) +{ + qmu=_qmu; +} NAMESPACE_END(Grid); diff --git a/tests/qdpxx/Test_qdpxx_munprec.cc b/tests/qdpxx/Test_qdpxx_munprec.cc index 82874546..c6ce2800 100644 --- a/tests/qdpxx/Test_qdpxx_munprec.cc +++ b/tests/qdpxx/Test_qdpxx_munprec.cc @@ -1,7 +1,6 @@ /************************************************************************************* Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc Copyright (C) 2015 @@ -26,13 +25,17 @@ Author: paboyle See the full license in the file "LICENSE" in the top level distribution directory *************************************************************************************/ /* END LEGAL */ +#include +#include +#include + #include int Ls=8; double M5=1.6; double mq=0.01; -double zolo_lo = 0.1; -double zolo_hi = 2.0; +double zolo_lo = 0.01; +double zolo_hi = 7.0; double mobius_scale=2.0; enum ChromaAction { @@ -55,11 +58,6 @@ enum ChromaAction { void calc_grid (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag); void calc_chroma (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag); -#include -#include -#include - - namespace Chroma { @@ -81,7 +79,7 @@ public: std::vector x(4); QDP::multi1d cx(4); - std::vector gd= gr.Grid()->GlobalDimensions(); + Grid::Coordinate gd = gr.Grid()->GlobalDimensions(); for (x[0]=0;x[0] x(5); QDP::multi1d cx(4); - std::vector gd= gr.Grid()->GlobalDimensions(); + Grid::Coordinate gd= gr.Grid()->GlobalDimensions(); for (x[0]=0;x[0] x(5); QDP::multi1d cx(4); - std::vector gd= gr.Grid()->GlobalDimensions(); + Grid::Coordinate gd= gr.Grid()->GlobalDimensions(); for (x[0]=0;x[0]OVEXT_CONSTANT_STRATEGY\n"; +"OVEXT_CONSTANT_STRATEGY1.0\n"; + UnprecOvExtFermActArray S_f(cfs,param); + Handle< FermState > fs( S_f.createState(u) ); + Handle< LinearOperatorArray > M(S_f.linOp(fs)); + return M; + } + if ( parms == HwPartFracTanh ) { + if ( Ls%2 == 0 ) { + printf("Ls is not odd\n"); + exit(-1); + } + UnprecOvExtFermActArrayParams param; + param.OverMass=M5; + param.Mass=_mq; + param.RatPolyDeg = Ls; + param.ApproxMin =eps_lo; + param.ApproxMax =eps_hi; + param.b5 =1.0; + param.c5 =1.0; + // param.approximation_type=COEFF_TYPE_ZOLOTAREV; + param.approximation_type=COEFF_TYPE_TANH_UNSCALED; + //param.approximation_type=COEFF_TYPE_TANH; + param.tuning_strategy_xml= + "OVEXT_CONSTANT_STRATEGY1.0\n"; UnprecOvExtFermActArray S_f(cfs,param); Handle< FermState > fs( S_f.createState(u) ); Handle< LinearOperatorArray > M(S_f.linOp(fs)); @@ -316,7 +337,35 @@ public: param.ApproxMin=eps_lo; param.ApproxMax=eps_hi; param.approximation_type=COEFF_TYPE_ZOLOTAREV; - param.RatPolyDeg=Ls; + param.RatPolyDeg=Ls-1; + // The following is why I think Chroma made some directional errors: + param.AuxFermAct= std::string( +"\n" +" UNPRECONDITIONED_WILSON\n" +" -1.8\n" +" 1\n" +" 0\n" +" 1000\n" +" 1.0e-9\n" +" \n" +" SIMPLE_FERMBC\n" +" 1 1 1 1\n" +" \n" +"" +); + param.AuxFermActGrp= std::string(""); + UnprecOvlapContFrac5DFermActArray S_f(fbc,param); + Handle< FermState > fs( S_f.createState(u) ); + Handle< LinearOperatorArray > M(S_f.linOp(fs)); + return M; + } + if ( parms == HwContFracTanh ) { + UnprecOvlapContFrac5DFermActParams param; + param.Mass=_mq; // How is M5 set? Wilson mass In AuxFermAct + param.ApproxMin=eps_lo; + param.ApproxMax=eps_hi; + param.approximation_type=COEFF_TYPE_TANH_UNSCALED; + param.RatPolyDeg=Ls-1; // The following is why I think Chroma made some directional errors: param.AuxFermAct= std::string( "\n" @@ -378,7 +427,14 @@ int main (int argc,char **argv ) * Setup QDP *********************************************************/ Chroma::initialize(&argc,&argv); - Chroma::WilsonTypeFermActs4DEnv::registerAll(); + // Chroma::WilsonTypeFermActs4DEnv::registerAll(); + Chroma::WilsonTypeFermActsEnv::registerAll(); + //bool linkageHack(void) + //{ + // bool foo = true; + // Inline Measurements + // InlineAggregateEnv::registerAll(); + // GaugeInitEnv::registerAll(); /******************************************************** * Setup Grid @@ -388,26 +444,34 @@ int main (int argc,char **argv ) Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()), Grid::GridDefaultMpi()); - std::vector gd = UGrid->GlobalDimensions(); + Grid::Coordinate gd = UGrid->GlobalDimensions(); QDP::multi1d nrow(QDP::Nd); for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu]; QDP::Layout::setLattSize(nrow); QDP::Layout::create(); - Grid::GridCartesian * FGrid = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); - Grid::LatticeGaugeField lat(UGrid); - Grid::LatticeFermion src(FGrid); - Grid::LatticeFermion res_chroma(FGrid); - Grid::LatticeFermion res_grid (FGrid); - std::vector ActionList({ HtCayleyTanh, // Plain old DWF. HmCayleyTanh, HwCayleyTanh, HtCayleyZolo, // Plain old DWF. HmCayleyZolo, - HwCayleyZolo + HwCayleyZolo, + HwPartFracZolo, + HwContFracZolo, + HwContFracTanh + }); + std::vector LsList({ + 8,//HtCayleyTanh, // Plain old DWF. + 8,//HmCayleyTanh, + 8,//HwCayleyTanh, + 8,//HtCayleyZolo, // Plain old DWF. + 8,//HmCayleyZolo, + 8,//HwCayleyZolo, + 9,//HwPartFracZolo + 9, //HwContFracZolo + 9 //HwContFracTanh }); std::vector ActionName({ "HtCayleyTanh", @@ -415,10 +479,19 @@ int main (int argc,char **argv ) "HwCayleyTanh", "HtCayleyZolo", "HmCayleyZolo", - "HwCayleyZolo" + "HwCayleyZolo", + "HwPartFracZolo", + "HwContFracZolo", + "HwContFracTanh" }); for(int i=0;i::HotConfiguration(RNG4,Umu); + Grid::SU::HotConfiguration(RNG4,Umu); /* Grid::LatticeColourMatrix U(UGrid); @@ -519,7 +593,7 @@ void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeF if ( action == HtCayleyTanh ) { - Grid::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5); + Grid::DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5); std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<