diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index be528e79..e47ff331 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -229,7 +229,14 @@ namespace QCD { } } - void CayleyFermion5D::SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c) + // Tanh + void CayleyFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c) + { + SetCoefficientsZolotarev(1.0,zdata,b,c); + + } + //Zolo + void CayleyFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c) { /////////////////////////////////////////////////////////// @@ -266,7 +273,7 @@ namespace QCD { double bmc = b-c; for(int i=0; i < Ls; i++){ as[i] = 1.0; - omega[i] = ((double)zdata->gamma[i]); //NB reciprocal relative to Chroma NEF code + omega[i] = ((double)zdata->gamma[i])*zolo_hi; //NB reciprocal relative to Chroma NEF code bs[i] = 0.5*(bpc/omega[i] + bmc); cs[i] = 0.5*(bpc/omega[i] - bmc); } diff --git a/lib/qcd/action/fermion/CayleyFermion5D.h b/lib/qcd/action/fermion/CayleyFermion5D.h index 57c71992..e2175d77 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.h +++ b/lib/qcd/action/fermion/CayleyFermion5D.h @@ -20,7 +20,7 @@ namespace Grid { virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); - + virtual void Instantiatable(void)=0; // protected: RealD mass; @@ -52,7 +52,8 @@ namespace Grid { RealD _mass,RealD _M5); protected: - void SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c); + void SetCoefficientsZolotarev(RealD zolohi,Approx::zolotarev_data *zdata,RealD b,RealD c); + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c); }; } diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc index c281b486..250e365f 100644 --- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.cc @@ -1,9 +1,56 @@ #include namespace Grid { - namespace QCD { + void ContinuedFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c) + { + SetCoefficientsZolotarev(1.0,zdata,b,c); + } + void ContinuedFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c) + { + R=(1+this->mass)/(1-this->mass); + + Beta.resize(Ls); + cc.resize(Ls); + cc_d.resize(Ls); + sqrt_cc.resize(Ls); + for(int i=0; i < Ls ; i++){ + Beta[i] = zdata -> beta[i]; + cc[i] = 1.0/Beta[i]; + cc_d[i]=sqrt(cc[i]); + } + + cc_d[Ls-1]=1.0; + for(int i=0; i < Ls-1 ; i++){ + sqrt_cc[i]= sqrt(cc[i]*cc[i+1]); + } + sqrt_cc[Ls-2]=sqrt(cc[Ls-2]); + + + ZoloHiInv =1.0/zolo_hi; + double dw_diag = (4.0-M5)*ZoloHiInv; + + See.resize(Ls); + Aee.resize(Ls); + int sign=1; + for(int s=0;sM5)*scale; + double dw_diag = (4.0-M5)*ZoloHiInv; int sign=1; for(int s=0;smass)/(1-this->mass); + double R=(1+mass)/(1-mass); ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1); ag5xpby_ssp(chi,R,psi,1.0,chi,s,s); } else { @@ -80,7 +131,7 @@ namespace Grid { void ContinuedFractionFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi) { // Apply Linv - axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); + axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); for(int s=1;sLs-1,this->Ls-1); + axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,Ls-1,Ls-1); for(int s=Ls-2;s>=0;s--){ axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1); } @@ -112,6 +163,10 @@ namespace Grid { FourDimGrid, FourDimRedBlackGrid,M5), mass(_mass) { + assert((Ls&0x1)==1); // Odd Ls required + int nrational=Ls-1;// Even rational order + zdata = Approx::grid_higham(1.0,nrational);// eps is ignored for higham + SetCoefficientsTanh(zdata,1.0,0.0); } } diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h index 7f5c022a..99365009 100644 --- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h @@ -21,20 +21,8 @@ namespace Grid { virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); - private: - - Approx::zolotarev_data *zdata; - - // Cont frac - RealD mass; - RealD R; - RealD scale; - std::vector Beta; - std::vector cc;; - std::vector cc_d;; - std::vector sqrt_cc; - std::vector See; - std::vector Aee; + // virtual void Instantiatable(void)=0; + virtual void Instantiatable(void) {}; // Constructors ContinuedFractionFermion5D(LatticeGaugeField &_Umu, @@ -44,6 +32,24 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD M5); + protected: + + void SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD b,RealD c); + void SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata,RealD b,RealD c); + + Approx::zolotarev_data *zdata; + + // Cont frac + RealD mass; + RealD R; + RealD ZoloHiInv; + std::vector Beta; + std::vector cc;; + std::vector cc_d;; + std::vector sqrt_cc; + std::vector See; + std::vector Aee; + }; diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index 3e6a9739..a25c0c3c 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -11,6 +11,7 @@ namespace Grid { { public: + virtual void Instantiatable(void) {}; // Constructors DomainWallFermion(LatticeGaugeField &_Umu, GridCartesian &FiveDimGrid, @@ -33,7 +34,7 @@ namespace Grid { std::cout << "DomainWallFermion with Ls="<CayleyFermion5D::SetCoefficients(1.0,zdata,1.0,0.0); + this->CayleyFermion5D::SetCoefficientsTanh(zdata,1.0,0.0); } diff --git a/lib/qcd/action/fermion/MobiusFermion.h b/lib/qcd/action/fermion/MobiusFermion.h index 4c291fad..33f94089 100644 --- a/lib/qcd/action/fermion/MobiusFermion.h +++ b/lib/qcd/action/fermion/MobiusFermion.h @@ -11,6 +11,7 @@ namespace Grid { { public: + virtual void Instantiatable(void) {}; // Constructors MobiusFermion(LatticeGaugeField &_Umu, GridCartesian &FiveDimGrid, @@ -34,7 +35,7 @@ namespace Grid { assert(zdata->n==this->Ls); // Call base setter - this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c); + this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c); } diff --git a/lib/qcd/action/fermion/MobiusZolotarevFermion.h b/lib/qcd/action/fermion/MobiusZolotarevFermion.h index 9ac795d9..1be61601 100644 --- a/lib/qcd/action/fermion/MobiusZolotarevFermion.h +++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h @@ -11,6 +11,7 @@ namespace Grid { { public: + virtual void Instantiatable(void) {}; // Constructors MobiusZolotarevFermion(LatticeGaugeField &_Umu, GridCartesian &FiveDimGrid, @@ -34,10 +35,9 @@ namespace Grid { assert(zdata->n==this->Ls); std::cout << "MobiusZolotarevFermion (b="<CayleyFermion5D::SetCoefficients(1.0,zdata,b,c); + this->CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c); } diff --git a/tests/InvSqrt.gnu b/tests/InvSqrt.gnu deleted file mode 100644 index e69de29b..00000000 diff --git a/tests/Make.inc b/tests/Make.inc index d592f218..b525874d 100644 --- a/tests/Make.inc +++ b/tests/Make.inc @@ -1,5 +1,21 @@ -bin_PROGRAMS = Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_many_cg Test_many_evenodd Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_evenodd +bin_PROGRAMS = Test_cayley_cg Test_cayley_even_odd Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_even_odd + + +Test_cayley_cg_SOURCES=Test_cayley_cg.cc +Test_cayley_cg_LDADD=-lGrid + + +Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc +Test_cayley_even_odd_LDADD=-lGrid + + +Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc +Test_contfrac_cg_LDADD=-lGrid + + +Test_contfrac_even_odd_SOURCES=Test_contfrac_even_odd.cc +Test_contfrac_even_odd_LDADD=-lGrid Test_cshift_SOURCES=Test_cshift.cc @@ -34,14 +50,6 @@ Test_main_SOURCES=Test_main.cc Test_main_LDADD=-lGrid -Test_many_cg_SOURCES=Test_many_cg.cc -Test_many_cg_LDADD=-lGrid - - -Test_many_evenodd_SOURCES=Test_many_evenodd.cc -Test_many_evenodd_LDADD=-lGrid - - Test_nersc_io_SOURCES=Test_nersc_io.cc Test_nersc_io_LDADD=-lGrid @@ -78,6 +86,6 @@ Test_wilson_cg_unprec_SOURCES=Test_wilson_cg_unprec.cc Test_wilson_cg_unprec_LDADD=-lGrid -Test_wilson_evenodd_SOURCES=Test_wilson_evenodd.cc -Test_wilson_evenodd_LDADD=-lGrid +Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc +Test_wilson_even_odd_LDADD=-lGrid diff --git a/tests/Sqrt.gnu b/tests/Sqrt.gnu deleted file mode 100644 index ae56ab97..00000000 --- a/tests/Sqrt.gnu +++ /dev/null @@ -1,2 +0,0 @@ -f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866)); -f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422)); diff --git a/tests/Test_many_cg.cc b/tests/Test_cayley_cg.cc similarity index 100% rename from tests/Test_many_cg.cc rename to tests/Test_cayley_cg.cc diff --git a/tests/Test_many_evenodd.cc b/tests/Test_cayley_even_odd.cc similarity index 100% rename from tests/Test_many_evenodd.cc rename to tests/Test_cayley_even_odd.cc diff --git a/tests/Test_contfrac_cg.cc b/tests/Test_contfrac_cg.cc new file mode 100644 index 00000000..7fa0d6fc --- /dev/null +++ b/tests/Test_contfrac_cg.cc @@ -0,0 +1,147 @@ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + std::cout <<"ContinuedFractionFermion test"<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} +template +void TestCGinversions(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + std::cout << "Testing unpreconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black preconditioned inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); + std::cout << "Testing red black Schur inverter"<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,RNG4,RNG5); +} + +template +void TestCGunprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + HermitianOperator HermOp(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOp,src,result); + +} +template +void TestCGprec(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion src_o(FrbGrid); + LatticeFermion result_o(FrbGrid); + pickCheckerboard(Odd,src_o,src); + result_o=zero; + + HermitianCheckerBoardedOperator HermOpEO(Ddwf); + ConjugateGradient CG(1.0e-8,10000); + CG(HermOpEO,src_o,result_o); +} + + +template +void TestCGschur(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion result(FGrid); result=zero; + + ConjugateGradient CG(1.0e-8,10000); + SchurRedBlackSolve SchurSolver(CG); + SchurSolver(Ddwf,src,result); +} diff --git a/tests/Test_contfrac_even_odd.cc b/tests/Test_contfrac_even_odd.cc new file mode 100644 index 00000000..801bd955 --- /dev/null +++ b/tests/Test_contfrac_even_odd.cc @@ -0,0 +1,218 @@ +#include + +using namespace std; +using namespace Grid; +using namespace Grid::QCD; + +template +struct scal { + d internal; +}; + + Gamma::GammaMatrix Gmu [] = { + Gamma::GammaX, + Gamma::GammaY, + Gamma::GammaZ, + Gamma::GammaT + }; + + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, GridParallelRNG *RNG5); + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + + int threads = GridThread::GetThreads(); + std::cout << "Grid is setup to use "< seeds4({1,2,3,4}); + std::vector seeds5({5,6,7,8}); + GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4); + + LatticeGaugeField Umu(UGrid); random(RNG4,Umu); + std::vector U(4,UGrid); + + RealD mass=0.1; + RealD M5 =1.8; + std::cout <<"ContinuedFractionFermion test"<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + Grid_finalize(); +} + +template +void TestWhat(What & Ddwf, + GridCartesian * FGrid, GridRedBlackCartesian * FrbGrid, + GridCartesian * UGrid, GridRedBlackCartesian * UrbGrid, + RealD mass, RealD M5, + GridParallelRNG *RNG4, + GridParallelRNG *RNG5) +{ + + LatticeFermion src (FGrid); random(*RNG5,src); + LatticeFermion phi (FGrid); random(*RNG5,phi); + LatticeFermion chi (FGrid); random(*RNG5,chi); + LatticeFermion result(FGrid); result=zero; + LatticeFermion ref(FGrid); ref=zero; + LatticeFermion tmp(FGrid); tmp=zero; + LatticeFermion err(FGrid); tmp=zero; + + LatticeFermion src_e (FrbGrid); + LatticeFermion src_o (FrbGrid); + LatticeFermion r_e (FrbGrid); + LatticeFermion r_o (FrbGrid); + LatticeFermion r_eo (FGrid); + LatticeFermion r_eeoo(FGrid); + + std::cout<<"=========================================================="< * = < chi | Deo^dag| phi> "<