From 59163862420a81a5e9d41bd11eae2c5c329e96b9 Mon Sep 17 00:00:00 2001 From: Peter Boyle Date: Wed, 3 Jun 2015 09:36:26 +0100 Subject: [PATCH] Mobius Caley form, Mobius Zolotarev operators. Pass Even Odd vs unprec test and hermiticity checks in tests/Grid_any_evenodd.cc; will work on inversion tests shortly. --- lib/qcd/action/Actions.h | 8 +- lib/qcd/action/fermion/CayleyFermion5D.cc | 110 ++++++++- lib/qcd/action/fermion/CayleyFermion5D.h | 6 +- lib/qcd/action/fermion/DomainWallFermion.h | 83 +------ lib/qcd/action/fermion/MobiusFermion.h | 46 ++++ .../action/fermion/MobiusZolotarevFermion.h | 48 ++++ lib/qcd/action/fermion/ScaledShamir.h | 51 ++++ lib/qcd/action/fermion/WilsonFermion5D.h | 12 + tests/Grid_any_evenodd.cc | 226 ++++++++++++++++++ tests/Makefile.am | 12 +- 10 files changed, 507 insertions(+), 95 deletions(-) create mode 100644 lib/qcd/action/fermion/MobiusFermion.h create mode 100644 lib/qcd/action/fermion/MobiusZolotarevFermion.h create mode 100644 lib/qcd/action/fermion/ScaledShamir.h create mode 100644 tests/Grid_any_evenodd.cc diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index acbf027c..d37d1cd4 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -38,6 +38,9 @@ //#include #include +#include +#include +#include //#include @@ -70,11 +73,6 @@ class LinearGaugeAction : public GaugeAction< multi1d, multi1d > typedef multi1d P; - typedef multi1d Q; - virtual void staple(LatticeColorMatrix& result, - const Handle< GaugeState >& state, - int mu, int cb) const = 0; */ - #endif diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index 263cc28b..be528e79 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -15,7 +15,6 @@ namespace QCD { FourDimRedBlackGrid,_M5), mass(_mass) { - std::cout << "Constructing a CayleyFermion5D"<gamma[i]); //NB reciprocal relative to Chroma NEF code + bs[i] = 0.5*(bpc/omega[i] + bmc); + cs[i] = 0.5*(bpc/omega[i] - bmc); + } + + //////////////////////////////////////////////////////// + // Constants for the preconditioned matrix Cayley form + //////////////////////////////////////////////////////// + bee.resize(Ls); + cee.resize(Ls); + beo.resize(Ls); + ceo.resize(Ls); + + for(int i=0;i omega; std::vector bs; // S dependent coeffs @@ -53,6 +51,8 @@ namespace Grid { GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5); + protected: + void SetCoefficients(RealD scale,Approx::zolotarev_data *zdata,RealD b,RealD c); }; } diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index 2abb6eb2..3e6a9739 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -28,86 +28,13 @@ namespace Grid { { RealD eps = 1.0; - zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham + Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham assert(zdata->n==this->Ls); + + std::cout << "DomainWallFermion with Ls="<CayleyFermion5D::SetCoefficients(1.0,zdata,1.0,0.0); - /////////////////////////////////////////////////////////// - // The Cayley coeffs (unprec) - /////////////////////////////////////////////////////////// - this->omega.resize(this->Ls); - this->bs.resize(this->Ls); - this->cs.resize(this->Ls); - this->as.resize(this->Ls); - - for(int i=0; i < this->Ls; i++){ - this->as[i] = 1.0; - this->omega[i] = ((double)zdata -> gamma[i]); - double bb=1.0; - this->bs[i] = 0.5*(bb/(this->omega[i]) + 1.0); - this->cs[i] = 0.5*(bb/(this->omega[i]) - 1.0); - } - - //////////////////////////////////////////////////////// - // Constants for the preconditioned matrix Cayley form - //////////////////////////////////////////////////////// - this->bee.resize(this->Ls); - this->cee.resize(this->Ls); - this->beo.resize(this->Ls); - this->ceo.resize(this->Ls); - - for(int i=0;iLs;i++){ - this->bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0); - this->cee[i]=as[i]*(1.0-cs[i]*(4.0-M5)); - this->beo[i]=as[i]*bs[i]; - this->ceo[i]=-as[i]*cs[i]; - } - - aee.resize(this->Ls); - aeo.resize(this->Ls); - for(int i=0;iLs;i++){ - aee[i]=cee[i]; - aeo[i]=ceo[i]; - } - - ////////////////////////////////////////// - // LDU decomposition of eeoo - ////////////////////////////////////////// - dee.resize(this->Ls); - lee.resize(this->Ls); - leem.resize(this->Ls); - uee.resize(this->Ls); - ueem.resize(this->Ls); - - for(int i=0;iLs;i++){ - - dee[i] = bee[i]; - - if ( i < this->Ls-1 ) { - - lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column - - leem[i]=this->mass*cee[this->Ls-1]/bee[0]; - for(int j=0;jmass; - for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j]; - ueem[i]*= aee[0]/bee[0]; - - } else { - lee[i] =0.0; - leem[i]=0.0; - uee[i] =0.0; - ueem[i]=0.0; - } - } - - { - double delta_d=mass*cee[this->Ls-1]; - for(int j=0;jLs-1;j++) delta_d *= cee[j]/bee[j]; - dee[this->Ls-1] += delta_d; - } } }; diff --git a/lib/qcd/action/fermion/MobiusFermion.h b/lib/qcd/action/fermion/MobiusFermion.h new file mode 100644 index 00000000..4c291fad --- /dev/null +++ b/lib/qcd/action/fermion/MobiusFermion.h @@ -0,0 +1,46 @@ +#ifndef GRID_QCD_MOBIUS_FERMION_H +#define GRID_QCD_MOBIUS_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class MobiusFermion : public CayleyFermion5D + { + public: + + // Constructors + MobiusFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD b, RealD c) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + + { + RealD eps = 1.0; + + std::cout << "MobiusFermion (b="<Ls);// eps is ignored for higham + assert(zdata->n==this->Ls); + + // Call base setter + this->CayleyFermion5D::SetCoefficients(1.0,zdata,b,c); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/MobiusZolotarevFermion.h b/lib/qcd/action/fermion/MobiusZolotarevFermion.h new file mode 100644 index 00000000..866d0c39 --- /dev/null +++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h @@ -0,0 +1,48 @@ +#ifndef GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H +#define GRID_QCD_MOBIUS_ZOLOTAREV_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class MobiusZolotarevFermion : public CayleyFermion5D + { + public: + + // Constructors + MobiusZolotarevFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, + RealD b, RealD c, + RealD lo, RealD hi) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) + + { + RealD eps = lo/hi; + + Approx::zolotarev_data *zdata = Approx::grid_zolotarev(eps,this->Ls,0);// eps is ignored for higham + assert(zdata->n==this->Ls); + + std::cout << "MobiusZolotarevFermion (b="<CayleyFermion5D::SetCoefficients(1.0,zdata,b,c); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/ScaledShamir.h b/lib/qcd/action/fermion/ScaledShamir.h new file mode 100644 index 00000000..a1fd33d0 --- /dev/null +++ b/lib/qcd/action/fermion/ScaledShamir.h @@ -0,0 +1,51 @@ +#ifndef GRID_QCD_DOMAIN_WALL_FERMION_H +#define GRID_QCD_DOMAIN_WALL_FERMION_H + +#include + +namespace Grid { + + namespace QCD { + + class ScaledShamirFermion : public CayleyFermion5D + { + public: + + // Constructors + ScaledShamirFermion(LatticeGaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5, RealD scale) : + + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5, + RealD b, + RealD c) + + { + RealD eps = 1.0; + + Approx::zolotarev_data *zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham + assert(zdata->n==this->Ls); + + //b+c = scale; + //b-c = 1 + //b = 0.5(scale+1); + //c = 0.5(scale-1); + + // Call base setter + this->CayleyFermion5D::SetCoefficients(1.0,zdata,0.5*(scale+1.0),0.5*(scale-1.0)); + + } + + }; + + } +} + +#endif diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d4777d01..062c4d82 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -14,6 +14,18 @@ namespace Grid { // i.e. even even contains fifth dim hopping term. // // [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ] + //////////////////////////// + //ContFrac: + // Ls always odd. Rational poly deg is either Ls or Ls-1 + //PartFrac + // Ls always odd. Rational poly deg is either Ls or Ls-1 + // + //Cayley: Ls always even, Rational poly deg is Ls + // + // Just set nrational as Ls. Forget about Ls-1 cases. + // + // Require odd Ls for cont and part frac + //////////////////////////// //////////////////////////////////////////////////////////////////////////////// class WilsonFermion5D : public FermionOperator { diff --git a/tests/Grid_any_evenodd.cc b/tests/Grid_any_evenodd.cc new file mode 100644 index 00000000..8d8580c1 --- /dev/null +++ b/tests/Grid_any_evenodd.cc @@ -0,0 +1,226 @@ +#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; + DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + TestWhat(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + RealD b=1.5;// Scale factor b+c=2, b-c=1 + RealD c=0.5; + MobiusFermion Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + TestWhat(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + + MobiusZolotarevFermion Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,4.0); + TestWhat(Dzolo,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> "< $f.log - done Grid_main_SOURCES = Grid_main.cc Grid_main_LDADD = -lGrid @@ -66,6 +63,9 @@ Grid_simd_LDADD = -lGrid Grid_wilson_evenodd_SOURCES = Grid_wilson_evenodd.cc Grid_wilson_evenodd_LDADD = -lGrid +Grid_any_evenodd_SOURCES = Grid_any_evenodd.cc +Grid_any_evenodd_LDADD = -lGrid + Grid_wilson_cg_unprec_SOURCES = Grid_wilson_cg_unprec.cc Grid_wilson_cg_unprec_LDADD = -lGrid