diff --git a/Grid/qcd/action/fermion/CayleyFermion5D.h b/Grid/qcd/action/fermion/CayleyFermion5D.h index c8fbe5a8..ec80b692 100644 --- a/Grid/qcd/action/fermion/CayleyFermion5D.h +++ b/Grid/qcd/action/fermion/CayleyFermion5D.h @@ -124,6 +124,11 @@ public: RealD _b; RealD _c; + // possible boost + std::vector qmu; + void set_qmu(std::vector _qmu) { qmu=_qmu; assert(qmu.size()==Nd);}; + void addQmu(const FermionField &in, FermionField &out, int dag); + // Cayley form Moebius (tanh and zolotarev) std::vector omega; std::vector bs; // S dependent coeffs diff --git a/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h b/Grid/qcd/action/fermion/ContinuedFractionFermion5D.h index bfc0dd8b..3fb84cd5 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 47406730..a71fc3f3 100644 --- a/Grid/qcd/action/fermion/PartialFractionFermion5D.h +++ b/Grid/qcd/action/fermion/PartialFractionFermion5D.h @@ -83,11 +83,70 @@ 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); + if ( this->qmu.size() ){ + this->MomentumSpacePropagatorHwQ(prop_k,in_k,mass,twist,this->qmu); + } else { + 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) { qmu=_qmu; assert(qmu.size()==Nd);}; + void addQmu(const FermionField &in, FermionField &out, int dag); + protected: virtual void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale); virtual void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata); + std::vector qmu; + // Part frac RealD mass; RealD dw_diag; diff --git a/Grid/qcd/action/fermion/WilsonFermion5D.h b/Grid/qcd/action/fermion/WilsonFermion5D.h index d614b290..dd83f269 100644 --- a/Grid/qcd/action/fermion/WilsonFermion5D.h +++ b/Grid/qcd/action/fermion/WilsonFermion5D.h @@ -109,6 +109,8 @@ public: void MomentumSpacePropagatorHt_5d(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; void MomentumSpacePropagatorHt(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; void MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass,std::vector twist) ; + void MomentumSpacePropagatorHwQ(FermionField &out,const FermionField &in,RealD mass,std::vector twist, + std::vector qmu) ; // Implement hopping term non-hermitian hopping term; half cb or both // Implement s-diagonal DW diff --git a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h index 69b5b02c..2ace6c18 100644 --- a/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/CayleyFermion5DImplementation.h @@ -48,7 +48,8 @@ CayleyFermion5D::CayleyFermion5D(GaugeField &_Umu, FourDimGrid, FourDimRedBlackGrid,_M5,p), mass_plus(_mass), mass_minus(_mass) -{ +{ + // qmu defaults to zero size; } /////////////////////////////////////////////////////////////// @@ -270,6 +271,34 @@ void CayleyFermion5D::MeooeDag5D (const FermionField &psi, FermionField M5Ddag(psi,psi,Din,lower,diag,upper); } +template +void CayleyFermion5D::addQmu(const FermionField &psi,FermionField &chi, int dag) +{ + if ( qmu.size() ) { + + Gamma::Algebra Gmu [] = { + Gamma::Algebra::GammaX, + Gamma::Algebra::GammaY, + Gamma::Algebra::GammaZ, + Gamma::Algebra::GammaT + }; + std::vector coeff(Nd); + ComplexD ci(0,1); + + assert(qmu.size()==Nd); + + for(int mu=0;mu void CayleyFermion5D::M (const FermionField &psi, FermionField &chi) { @@ -277,8 +306,12 @@ void CayleyFermion5D::M (const FermionField &psi, FermionField &chi) // Assemble Din Meooe5D(psi,Din); - + this->DW(Din,chi,DaggerNo); + + // add i q_mu gamma_mu here + addQmu(Din,chi,DaggerNo); + // ((b D_W + D_w hop terms +1) on s-diag axpby(chi,1.0,1.0,chi,psi); @@ -295,6 +328,9 @@ void CayleyFermion5D::Mdag (const FermionField &psi, FermionField &chi) FermionField Din(psi.Grid()); // Apply Dw this->DW(psi,Din,DaggerYes); + + // add -i conj(q_mu) gamma_mu here ... if qmu is real, gammm_5 hermitian, otherwise not. + addQmu(psi,Din,DaggerYes); MeooeDag5D(Din,chi); 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 93684929..84884c6d 100644 --- a/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/PartialFractionFermion5DImplementation.h @@ -237,7 +237,32 @@ void PartialFractionFermion5D::M_internal(const FermionField &psi, Fermi // ( 0 -sqrt(p_i)*amax | 2 R gamma_5 + p0/amax 2H // - this->DW(psi,D,DaggerNo); + this->DW(psi,D,DaggerNo); + + // DW - DW+iqslash + // (g5 Dw)^dag = g5 Dw + // (iqmu g5 gmu)^dag = (-i qmu gmu^dag g5^dag) = i qmu g5 gmu + if ( qmu.size() ) { + + std::cout<< "Mat" << "qmu ("<::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 + //R g5 psi[Ls] + p[0] Hw ag5xpbg5y_ssp(chi,R*scale,psi,p[nblock]*scale/amax,D,Ls-1,Ls-1); - + for(int b=0;b::PartialFractionFermion5D(GaugeField &_Umu, { int Ls = this->Ls; - + qmu.resize(0); assert((Ls&0x1)==1); // Odd Ls required int nrational=Ls-1; @@ -461,6 +526,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/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h b/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h index 32e8108e..ec3bd94a 100644 --- a/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h +++ b/Grid/qcd/action/fermion/implementation/WilsonFermion5DImplementation.h @@ -740,6 +740,15 @@ void WilsonFermion5D::MomentumSpacePropagatorHt(FermionField &out,const Fe template void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const FermionField &in,RealD mass,std::vector twist) +{ + std::vector empty_q(Nd,0.0); + MomentumSpacePropagatorHwQ(out,in,mass,twist,empty_q); +} +template +void WilsonFermion5D::MomentumSpacePropagatorHwQ(FermionField &out,const FermionField &in, + RealD mass, + std::vector twist, + std::vector qmu) { Gamma::Algebra Gmu [] = { Gamma::Algebra::GammaX, @@ -755,6 +764,7 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe typedef typename FermionField::scalar_type ScalComplex; typedef Lattice > LatComplex; + typedef iSpinMatrix SpinMat; Coordinate latt_size = _grid->_fdimensions; @@ -772,8 +782,10 @@ void WilsonFermion5D::MomentumSpacePropagatorHw(FermionField &out,const Fe LatComplex kmu(_grid); ScalComplex ci(0.0,1.0); + std::cout<< "Feynman Rule" << "qmu ("<::MomentumSpacePropagatorHw(FermionField &out,const Fe kmu = kmu + TwoPiL * one * twist[mu];//momentum for twisted boundary conditions sk2 = sk2 + 2.0*sin(kmu*0.5)*sin(kmu*0.5); - sk = sk + sin(kmu)*sin(kmu); - num = num - sin(kmu)*ci*(Gamma(Gmu[mu])*in); + sk = sk + (sin(kmu)+qmu[mu])*(sin(kmu)+qmu[mu]); + + // Terms for boosted Fermion + // 1/2 [ -i gamma.(sin p + q ) ] + // [ --------------------- + 1 ] + // [ wq + b ] + // + // wq = sqrt( (sinp+q)^2 + b^2 ) + // + + num = num - (sin(kmu)+qmu[mu])*ci*(Gamma(Gmu[mu])*in); } num = num + mass * in ; diff --git a/tests/core/Test_fft_pf.cc b/tests/core/Test_fft_pf.cc index d60f37ee..52540e57 100644 --- a/tests/core/Test_fft_pf.cc +++ b/tests/core/Test_fft_pf.cc @@ -74,7 +74,7 @@ int main (int argc, char ** argv) { std::cout<<"****************************************"< HermOp(Dov); + MdagMLinearOperator HermOp(Dov); ConjugateGradient CG(1.0e-8,10000); CG(HermOp,src5,result5); std::cout << " Solved by Conjugate Gradient (CGNE)" < qmu({1.0,0.0,0.0,0.0}); + Dov.set_qmu(qmu); // Momentum space prop std::cout << " Solving by FFT and Feynman rules" < HermOp(Dov); + MdagMLinearOperator HermOp(Dov); ConjugateGradient CG(1.0e-8,10000); CG(HermOp,src5,result5); //////////////////////////////////////////////////////////////////////// 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 "<::HotConfiguration(RNG4, Umu); +#endif std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << " Ls: " << Ls << std::endl; diff --git a/tests/solver/Test_dwf_cg_unprec.cc b/tests/solver/Test_dwf_cg_unprec.cc index 58614c49..7435bfae 100644 --- a/tests/solver/Test_dwf_cg_unprec.cc +++ b/tests/solver/Test_dwf_cg_unprec.cc @@ -54,15 +54,30 @@ int main (int argc, char ** argv) GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid); GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid); + std::vector qmu; + qmu.push_back(ComplexD(0.1,0.0)); + qmu.push_back(ComplexD(0.0,0.0)); + qmu.push_back(ComplexD(0.0,0.0)); + qmu.push_back(ComplexD(0.0,0.01)); + + std::vector seeds4({1,2,3,4}); std::vector seeds5({5,6,7,8}); GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + LatticeFermion tmp(FGrid); LatticeFermion src(FGrid); random(RNG5,src); LatticeFermion result(FGrid); result=Zero(); - LatticeGaugeField Umu(UGrid); SU::HotConfiguration(RNG4,Umu); - + LatticeGaugeField Umu(UGrid); +#if 0 + FieldMetaData header; + std::string file("ckpoint_lat.4000"); + NerscIO::readConfiguration(Umu,header,file); +#else + SU::HotConfiguration(RNG4,Umu); +#endif + std::vector U(4,UGrid); for(int mu=0;mu(Umu,mu); @@ -71,8 +86,15 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.8; DomainWallFermionD Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + Ddwf.qmu = qmu; + Ddwf.M(src,tmp); + std::cout << " |M src|^2 "< HermOp(Ddwf); + HermOp.HermOp(src,tmp); + + std::cout << " "< CG(1.0e-6,10000); CG(HermOp,src,result);