diff --git a/benchmarks/Benchmark_dwf.cc b/benchmarks/Benchmark_dwf.cc index 9d9303ae..27772fdb 100644 --- a/benchmarks/Benchmark_dwf.cc +++ b/benchmarks/Benchmark_dwf.cc @@ -79,7 +79,7 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5 =1.8; - DomainWallFermion Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); std::cout<,SpinorIndex>::value,iVector >::type *SFINAE; + + const int SpinorIndex = 2; + template struct isSpinor { + static const bool value = (SpinorIndex==T::TensorLevel); + }; + template using IfSpinor = Invoke::value,int> > ; + template using IfNotSpinor = Invoke::value,int> > ; // ChrisK very keen to add extra space for Gparity doubling. // @@ -49,6 +59,10 @@ namespace QCD { template using iHalfSpinVector = iScalar, Nhs> >; template using iHalfSpinColourVector = iScalar, Nhs> >; + template using iGparitySpinColourVector = iVector, Nhs>, Ngp >; + template using iGparityHalfSpinColourVector = iVector, Nhs>, Ngp >; + + // Spin matrix typedef iSpinMatrix SpinMatrix; typedef iSpinMatrix SpinMatrixF; diff --git a/lib/qcd/action/Actions.h b/lib/qcd/action/Actions.h index f40cd554..9b119858 100644 --- a/lib/qcd/action/Actions.h +++ b/lib/qcd/action/Actions.h @@ -17,9 +17,6 @@ //////////////////////////////////////////// #include -#include - - //////////////////////////////////////////// // Gauge Actions //////////////////////////////////////////// @@ -32,58 +29,118 @@ #include //used by all wilson type fermions #include //used by all wilson type fermions + +//////////////////////////////////////////////////////////////////////////////////////////////////// +// Explicit explicit template instantiation is still required in the .cc files +// +// - CayleyFermion5D.cc +// - PartialFractionFermion5D.cc +// - WilsonFermion5D.cc +// - WilsonKernelsHand.cc +// - ContinuedFractionFermion5D.cc +// - WilsonFermion.cc +// - WilsonKernels.cc +// +// The explicit instantiation is only avoidable if we move this source to headers and end up with include/parse/recompile +// for EVERY .cc file. This define centralises the list and restores global push of impl cases +//////////////////////////////////////////////////////////////////////////////////////////////////// + +#define FermOpTemplateInstantiate(A) \ + template class A; \ + template class A; + //////////////////////////////////////////// -// 4D formulations +// Fermion operators / actions //////////////////////////////////////////// -#include +#include + +#include // 4d wilson like //#include -//////////////////////////////////////////// -// 5D formulations... -//////////////////////////////////////////// - -#include // used by all 5d overlap types - -////////// -// Cayley -////////// -#include +#include // 5d base used by all 5d overlap types +#include // Cayley types #include #include - #include #include -#include - #include #include +#include #include -////////////////////// -// Continued fraction -////////////////////// -#include -#include -#include +#include // Continued fraction +#include +#include -////////////////////// -// Partial fraction -////////////////////// -#include +#include // Partial fraction #include #include +//////////////////////////////////////////////////////////////////////////////////////////////////// +// More maintainable to maintain the following typedef list centrally, as more "impl" targets +// are added, (e.g. extension for gparity, half precision project in comms etc..) +//////////////////////////////////////////////////////////////////////////////////////////////////// + + +// Cayley 5d +namespace Grid { + namespace QCD { + +typedef WilsonFermion WilsonFermionR; +typedef WilsonFermion WilsonFermionF; +typedef WilsonFermion WilsonFermionD; + +typedef DomainWallFermion DomainWallFermionR; +typedef DomainWallFermion DomainWallFermionF; +typedef DomainWallFermion DomainWallFermionD; +typedef MobiusFermion MobiusFermionR; +typedef MobiusFermion MobiusFermionF; +typedef MobiusFermion MobiusFermionD; +typedef ScaledShamirFermion ScaledShamirFermionR; +typedef ScaledShamirFermion ScaledShamirFermionF; +typedef ScaledShamirFermion ScaledShamirFermionD; + +typedef MobiusZolotarevFermion MobiusZolotarevFermionR; +typedef MobiusZolotarevFermion MobiusZolotarevFermionF; +typedef MobiusZolotarevFermion MobiusZolotarevFermionD; +typedef ShamirZolotarevFermion ShamirZolotarevFermionR; +typedef ShamirZolotarevFermion ShamirZolotarevFermionF; +typedef ShamirZolotarevFermion ShamirZolotarevFermionD; + +typedef OverlapWilsonCayleyTanhFermion OverlapWilsonCayleyTanhFermionR; +typedef OverlapWilsonCayleyTanhFermion OverlapWilsonCayleyTanhFermionF; +typedef OverlapWilsonCayleyTanhFermion OverlapWilsonCayleyTanhFermionD; +typedef OverlapWilsonCayleyZolotarevFermion OverlapWilsonCayleyZolotarevFermionR; +typedef OverlapWilsonCayleyZolotarevFermion OverlapWilsonCayleyZolotarevFermionF; +typedef OverlapWilsonCayleyZolotarevFermion OverlapWilsonCayleyZolotarevFermionD; + +// Continued fraction +typedef OverlapWilsonContFracTanhFermion OverlapWilsonContFracTanhFermionR; +typedef OverlapWilsonContFracTanhFermion OverlapWilsonContFracTanhFermionF; +typedef OverlapWilsonContFracTanhFermion OverlapWilsonContFracTanhFermionD; +typedef OverlapWilsonContFracZolotarevFermion OverlapWilsonContFracZolotarevFermionR; +typedef OverlapWilsonContFracZolotarevFermion OverlapWilsonContFracZolotarevFermionF; +typedef OverlapWilsonContFracZolotarevFermion OverlapWilsonContFracZolotarevFermionD; + +// Partial fraction +typedef OverlapWilsonPartialFractionTanhFermion OverlapWilsonPartialFractionTanhFermionR; +typedef OverlapWilsonPartialFractionTanhFermion OverlapWilsonPartialFractionTanhFermionF; +typedef OverlapWilsonPartialFractionTanhFermion OverlapWilsonPartialFractionTanhFermionD; + +typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionR; +typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionF; +typedef OverlapWilsonPartialFractionZolotarevFermion OverlapWilsonPartialFractionZolotarevFermionD; + }} /////////////////////////////////////////////////////////////////////////////// // G5 herm -- this has to live in QCD since dirac matrix is not in the broader sector of code /////////////////////////////////////////////////////////////////////////////// #include //////////////////////////////////////// -// Pseudo fermion combinations +// Pseudo fermion combinations for HMC //////////////////////////////////////// #include #include - #endif diff --git a/lib/qcd/action/fermion/CayleyFermion5D.cc b/lib/qcd/action/fermion/CayleyFermion5D.cc index c18edb7b..1ec7f930 100644 --- a/lib/qcd/action/fermion/CayleyFermion5D.cc +++ b/lib/qcd/action/fermion/CayleyFermion5D.cc @@ -2,13 +2,14 @@ namespace Grid { namespace QCD { - CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD _M5) : - WilsonFermion5D(_Umu, + template + CayleyFermion5D::CayleyFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD _M5) : + WilsonFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, @@ -17,9 +18,11 @@ namespace QCD { { } - void CayleyFermion5D::Meooe5D (const LatticeFermion &psi, LatticeFermion &Din) + template + void CayleyFermion5D::Meooe5D (const FermionField &psi, FermionField &Din) { // Assemble Din + int Ls=this->Ls; for(int s=0;s + void CayleyFermion5D::MeooeDag5D (const FermionField &psi, FermionField &Din) { + int Ls=this->Ls; for(int s=0;s + RealD CayleyFermion5D::M (const FermionField &psi, FermionField &chi) { - LatticeFermion Din(psi._grid); + int Ls=this->Ls; + + FermionField Din(psi._grid); // Assemble Din /* @@ -75,7 +83,7 @@ namespace QCD { */ Meooe5D(psi,Din); - DW(Din,chi,DaggerNo); + this->DW(Din,chi,DaggerNo); // ((b D_W + D_w hop terms +1) on s-diag axpby(chi,1.0,1.0,chi,psi); @@ -95,18 +103,20 @@ namespace QCD { return norm2(chi); } - RealD CayleyFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi) + template + RealD CayleyFermion5D::Mdag (const FermionField &psi, FermionField &chi) { // Under adjoint //D1+ D1- P- -> D1+^dag P+ D2-^dag //D2- P+ D2+ P-D1-^dag D2+dag - LatticeFermion Din(psi._grid); + FermionField Din(psi._grid); // Apply Dw - DW(psi,Din,DaggerYes); + this->DW(psi,Din,DaggerYes); Meooe5D(Din,chi); + int Ls=this->Ls; for(int s=0;s + void CayleyFermion5D::Meooe (const FermionField &psi, FermionField &chi) { - LatticeFermion tmp(psi._grid); + FermionField tmp(psi._grid); // Assemble the 5d matrix Meooe5D(psi,tmp); std::cout << "Meooe Test replacement norm2 tmp = " <Ls; for(int s=0;sdb << " -db"<dn << " -dn"<dd << " -dd"<Ls; assert(zdata->db==Ls);// Beta has Ls coeffs R=(1+this->mass)/(1-this->mass); @@ -39,7 +41,7 @@ namespace Grid { ZoloHiInv =1.0/zolo_hi; - dw_diag = (4.0-M5)*ZoloHiInv; + dw_diag = (4.0-this->M5)*ZoloHiInv; See.resize(Ls); Aee.resize(Ls); @@ -61,11 +63,14 @@ namespace Grid { - RealD ContinuedFractionFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi) + template + RealD ContinuedFractionFermion5D::M (const FermionField &psi, FermionField &chi) { - LatticeFermion D(psi._grid); + int Ls = this->Ls; - DW(psi,D,DaggerNo); + FermionField D(psi._grid); + + this->DW(psi,D,DaggerNo); int sign=1; for(int s=0;s + RealD ContinuedFractionFermion5D::Mdag (const FermionField &psi, FermionField &chi) { // This matrix is already hermitian. (g5 Dw) = Dw dag g5 = (g5 Dw)dag // The rest of matrix is symmetric. // Can ignore "dag" return M(psi,chi); } - void ContinuedFractionFermion5D::Mdir (const LatticeFermion &psi, LatticeFermion &chi,int dir,int disp){ - DhopDir(psi,chi,dir,disp); // Dslash on diagonal. g5 Dslash is hermitian + template + void ContinuedFractionFermion5D::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){ + int Ls = this->Ls; + + this->DhopDir(psi,chi,dir,disp); // Dslash on diagonal. g5 Dslash is hermitian + int sign=1; for(int s=0;s + void ContinuedFractionFermion5D::Meooe (const FermionField &psi, FermionField &chi) { + int Ls = this->Ls; + // Apply 4d dslash if ( psi.checkerboard == Odd ) { - DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian + this->DhopEO(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian } else { - DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian + this->DhopOE(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian } int sign=1; @@ -121,12 +134,16 @@ namespace Grid { sign=-sign; } } - void ContinuedFractionFermion5D::MeooeDag (const LatticeFermion &psi, LatticeFermion &chi) + template + void ContinuedFractionFermion5D::MeooeDag (const FermionField &psi, FermionField &chi) { - Meooe(psi,chi); + this->Meooe(psi,chi); } - void ContinuedFractionFermion5D::Mooee (const LatticeFermion &psi, LatticeFermion &chi) + template + void ContinuedFractionFermion5D::Mooee (const FermionField &psi, FermionField &chi) { + int Ls = this->Ls; + int sign=1; for(int s=0;s + void ContinuedFractionFermion5D::MooeeDag (const FermionField &psi, FermionField &chi) { - Mooee(psi,chi); + this->Mooee(psi,chi); } - void ContinuedFractionFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi) + template + void ContinuedFractionFermion5D::MooeeInv (const FermionField &psi, FermionField &chi) { + int Ls = this->Ls; + // Apply Linv axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0); for(int s=1;s + void ContinuedFractionFermion5D::MooeeInvDag (const FermionField &psi, FermionField &chi) { - MooeeInv(psi,chi); + this->MooeeInv(psi,chi); } // force terms; five routines; default to Dhop on diagonal - void ContinuedFractionFermion5D::MDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void ContinuedFractionFermion5D::MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int sign=1; for(int s=0;sDhopDeriv(mat,D,V,DaggerNo); }; - void ContinuedFractionFermion5D::MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void ContinuedFractionFermion5D::MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int sign=1; for(int s=0;sDhopDerivOE(mat,D,V,DaggerNo); }; - void ContinuedFractionFermion5D::MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void ContinuedFractionFermion5D::MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int sign=1; for(int s=0;sDhopDerivEO(mat,D,V,DaggerNo); }; // Constructors - ContinuedFractionFermion5D::ContinuedFractionFermion5D( - LatticeGaugeField &_Umu, + template + ContinuedFractionFermion5D::ContinuedFractionFermion5D( + GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD M5) : - WilsonFermion5D(_Umu, + WilsonFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,M5), mass(_mass) { + int Ls = this->Ls; assert((Ls&0x1)==1); // Odd Ls required } + FermOpTemplateInstantiate(ContinuedFractionFermion5D); + } } diff --git a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h index 7d0c75bc..b1a7fba2 100644 --- a/lib/qcd/action/fermion/ContinuedFractionFermion5D.h +++ b/lib/qcd/action/fermion/ContinuedFractionFermion5D.h @@ -5,35 +5,38 @@ namespace Grid { namespace QCD { - class ContinuedFractionFermion5D : public WilsonFermion5D + template + class ContinuedFractionFermion5D : public WilsonFermion5D { public: +#include + public: // override multiply - virtual RealD M (const LatticeFermion &in, LatticeFermion &out); - virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); + virtual RealD M (const FermionField &in, FermionField &out); + virtual RealD Mdag (const FermionField &in, FermionField &out); // half checkerboard operaions - virtual void Meooe (const LatticeFermion &in, LatticeFermion &out); - virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out); - virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); - 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 Meooe (const FermionField &in, FermionField &out); + virtual void MeooeDag (const FermionField &in, FermionField &out); + virtual void Mooee (const FermionField &in, FermionField &out); + virtual void MooeeDag (const FermionField &in, FermionField &out); + virtual void MooeeInv (const FermionField &in, FermionField &out); + virtual void MooeeInvDag (const FermionField &in, FermionField &out); // force terms; five routines; default to Dhop on diagonal - virtual void MDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); // virtual void Instantiatable(void)=0; virtual void Instantiatable(void) =0; // Efficient support for multigrid coarsening - virtual void Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp); + virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp); // Constructors - ContinuedFractionFermion5D(LatticeGaugeField &_Umu, + ContinuedFractionFermion5D(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, diff --git a/lib/qcd/action/fermion/DomainWallFermion.h b/lib/qcd/action/fermion/DomainWallFermion.h index 498d8b34..00a6baba 100644 --- a/lib/qcd/action/fermion/DomainWallFermion.h +++ b/lib/qcd/action/fermion/DomainWallFermion.h @@ -7,24 +7,27 @@ namespace Grid { namespace QCD { - class DomainWallFermion : public CayleyFermion5D + template + class DomainWallFermion : public CayleyFermion5D { public: +#include + public: virtual void Instantiatable(void) {}; // Constructors - DomainWallFermion(LatticeGaugeField &_Umu, + DomainWallFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD _M5) : - CayleyFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + CayleyFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) { RealD eps = 1.0; @@ -32,9 +35,9 @@ namespace Grid { Approx::zolotarev_data *zdata = Approx::higham(eps,this->Ls);// eps is ignored for higham assert(zdata->n==this->Ls); - std::cout<Ls);// eps is ignored for higham assert(zdata->n==this->Ls); // Call base setter - this->CayleyFermion5D::SetCoefficientsTanh(zdata,b,c); + this->SetCoefficientsTanh(zdata,b,c); Approx::zolotarev_free(zdata); diff --git a/lib/qcd/action/fermion/MobiusZolotarevFermion.h b/lib/qcd/action/fermion/MobiusZolotarevFermion.h index 5f74702c..1c572326 100644 --- a/lib/qcd/action/fermion/MobiusZolotarevFermion.h +++ b/lib/qcd/action/fermion/MobiusZolotarevFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class MobiusZolotarevFermion : public CayleyFermion5D + template + class MobiusZolotarevFermion : public CayleyFermion5D { public: +#include + public: virtual void Instantiatable(void) {}; // Constructors - MobiusZolotarevFermion(LatticeGaugeField &_Umu, + MobiusZolotarevFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -22,7 +25,7 @@ namespace Grid { RealD b, RealD c, RealD lo, RealD hi) : - CayleyFermion5D(_Umu, + CayleyFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, @@ -34,10 +37,10 @@ namespace Grid { Approx::zolotarev_data *zdata = Approx::zolotarev(eps,this->Ls,0); assert(zdata->n==this->Ls); - std::cout<CayleyFermion5D::SetCoefficientsZolotarev(hi,zdata,b,c); + this->SetCoefficientsZolotarev(hi,zdata,b,c); Approx::zolotarev_free(zdata); } diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h index e764c8ae..0e8d4e19 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyTanhFermion.h @@ -7,12 +7,15 @@ namespace Grid { namespace QCD { - class OverlapWilsonCayleyTanhFermion : public MobiusFermion + template + class OverlapWilsonCayleyTanhFermion : public MobiusFermion { public: +#include + public: // Constructors - OverlapWilsonCayleyTanhFermion(LatticeGaugeField &_Umu, + OverlapWilsonCayleyTanhFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -21,11 +24,11 @@ namespace Grid { RealD scale) : // b+c=scale, b-c = 0 <=> b =c = scale/2 - MobiusFermion(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale) + MobiusFermion(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5,0.5*scale,0.5*scale) { } }; diff --git a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h index 82c43fb7..549292db 100644 --- a/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonCayleyZolotarevFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion + template + class OverlapWilsonCayleyZolotarevFermion : public MobiusZolotarevFermion { public: +#include + public: // Constructors - OverlapWilsonCayleyZolotarevFermion(LatticeGaugeField &_Umu, + OverlapWilsonCayleyZolotarevFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -21,7 +24,7 @@ namespace Grid { RealD _mass,RealD _M5, RealD lo, RealD hi) : // b+c=1.0, b-c = 0 <=> b =c = 1/2 - MobiusZolotarevFermion(_Umu, + MobiusZolotarevFermion(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h index 6ecd7822..73baffa7 100644 --- a/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonContfracTanhFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D + template + class OverlapWilsonContFracTanhFermion : public ContinuedFractionFermion5D { public: +#include + public: virtual void Instantiatable(void){}; // Constructors - OverlapWilsonContFracTanhFermion(LatticeGaugeField &_Umu, + OverlapWilsonContFracTanhFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -22,16 +25,16 @@ namespace Grid { RealD scale) : // b+c=scale, b-c = 0 <=> b =c = scale/2 - ContinuedFractionFermion5D(_Umu, + ContinuedFractionFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,_mass,_M5) { - assert((Ls&0x1)==1); // Odd Ls required - int nrational=Ls-1;// Even rational order + assert((this->Ls&0x1)==1); // Odd Ls required + int nrational=this->Ls-1;// Even rational order Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham - SetCoefficientsTanh(zdata,scale); + this->SetCoefficientsTanh(zdata,scale); Approx::zolotarev_free(zdata); } }; diff --git a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h index 595daef3..3c662e92 100644 --- a/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonContfracZolotarevFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D + template + class OverlapWilsonContFracZolotarevFermion : public ContinuedFractionFermion5D { public: +#include + public: virtual void Instantiatable(void){}; // Constructors - OverlapWilsonContFracZolotarevFermion(LatticeGaugeField &_Umu, + OverlapWilsonContFracZolotarevFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -22,19 +25,19 @@ namespace Grid { RealD lo,RealD hi): // b+c=scale, b-c = 0 <=> b =c = scale/2 - ContinuedFractionFermion5D(_Umu, + ContinuedFractionFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,_mass,_M5) { - assert((Ls&0x1)==1); // Odd Ls required + assert((this->Ls&0x1)==1); // Odd Ls required - int nrational=Ls;// Odd rational order + int nrational=this->Ls;// Odd rational order RealD eps = lo/hi; Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0); - SetCoefficientsZolotarev(hi,zdata); + this->SetCoefficientsZolotarev(hi,zdata); Approx::zolotarev_free(zdata); } diff --git a/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h b/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h index 1fa6c099..b017a9d7 100644 --- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionTanhFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D + template + class OverlapWilsonPartialFractionTanhFermion : public PartialFractionFermion5D { public: +#include + public: virtual void Instantiatable(void){}; // Constructors - OverlapWilsonPartialFractionTanhFermion(LatticeGaugeField &_Umu, + OverlapWilsonPartialFractionTanhFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -22,16 +25,16 @@ namespace Grid { RealD scale) : // b+c=scale, b-c = 0 <=> b =c = scale/2 - PartialFractionFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + PartialFractionFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) { - assert((Ls&0x1)==1); // Odd Ls required - int nrational=Ls-1;// Even rational order + assert((this->Ls&0x1)==1); // Odd Ls required + int nrational=this->Ls-1;// Even rational order Approx::zolotarev_data *zdata = Approx::higham(1.0,nrational);// eps is ignored for higham - SetCoefficientsTanh(zdata,scale); + this->SetCoefficientsTanh(zdata,scale); Approx::zolotarev_free(zdata); } }; diff --git a/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h b/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h index 3997c17f..3294f079 100644 --- a/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h +++ b/lib/qcd/action/fermion/OverlapWilsonPartialFractionZolotarevFermion.h @@ -7,13 +7,16 @@ namespace Grid { namespace QCD { - class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D + template + class OverlapWilsonPartialFractionZolotarevFermion : public PartialFractionFermion5D { public: +#include + public: virtual void Instantiatable(void){}; // Constructors - OverlapWilsonPartialFractionZolotarevFermion(LatticeGaugeField &_Umu, + OverlapWilsonPartialFractionZolotarevFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -22,19 +25,19 @@ namespace Grid { RealD lo,RealD hi): // b+c=scale, b-c = 0 <=> b =c = scale/2 - PartialFractionFermion5D(_Umu, - FiveDimGrid, - FiveDimRedBlackGrid, - FourDimGrid, - FourDimRedBlackGrid,_mass,_M5) + PartialFractionFermion5D(_Umu, + FiveDimGrid, + FiveDimRedBlackGrid, + FourDimGrid, + FourDimRedBlackGrid,_mass,_M5) { - assert((Ls&0x1)==1); // Odd Ls required + assert((this->Ls&0x1)==1); // Odd Ls required - int nrational=Ls;// Odd rational order + int nrational=this->Ls;// Odd rational order RealD eps = lo/hi; Approx::zolotarev_data *zdata = Approx::zolotarev(eps,nrational,0); - SetCoefficientsZolotarev(hi,zdata); + this->SetCoefficientsZolotarev(hi,zdata); Approx::zolotarev_free(zdata); } diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.cc b/lib/qcd/action/fermion/PartialFractionFermion5D.cc index 74174f51..81afb653 100644 --- a/lib/qcd/action/fermion/PartialFractionFermion5D.cc +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.cc @@ -2,12 +2,15 @@ namespace Grid { namespace QCD { - void PartialFractionFermion5D::Mdir (const LatticeFermion &psi, LatticeFermion &chi,int dir,int disp){ + + template + void PartialFractionFermion5D::Mdir (const FermionField &psi, FermionField &chi,int dir,int disp){ // this does both dag and undag but is trivial; make a common helper routing int sign = 1; + int Ls = this->Ls; - DhopDir(psi,chi,dir,disp); + this->DhopDir(psi,chi,dir,disp); int nblock=(Ls-1)/2; for(int b=0;b + void PartialFractionFermion5D::Meooe_internal(const FermionField &psi, FermionField &chi,int dag) { - // this does both dag and undag but is trivial; make a common helper routing + int Ls = this->Ls; int sign = dag ? (-1) : 1; if ( psi.checkerboard == Odd ) { - DhopEO(psi,chi,DaggerNo); + this->DhopEO(psi,chi,DaggerNo); } else { - DhopOE(psi,chi,DaggerNo); + this->DhopOE(psi,chi,DaggerNo); } int nblock=(Ls-1)/2; @@ -38,10 +42,12 @@ namespace Grid { ag5xpby_ssp(chi,p[nblock]*scale/amax,chi,0.0,chi,Ls-1,Ls-1); } - void PartialFractionFermion5D::Mooee_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag) + template + void PartialFractionFermion5D::Mooee_internal(const FermionField &psi, FermionField &chi,int dag) { // again dag and undag are trivially related int sign = dag ? (-1) : 1; + int Ls = this->Ls; int nblock=(Ls-1)/2; for(int b=0;b + void PartialFractionFermion5D::MooeeInv_internal(const FermionField &psi, FermionField &chi,int dag) { int sign = dag ? (-1) : 1; + int Ls = this->Ls; - LatticeFermion tmp(psi._grid); + FermionField tmp(psi._grid); /////////////////////////////////////////////////////////////////////////////////////// //Linv @@ -129,10 +137,12 @@ namespace Grid { axpby_ssp (chi, 1.0/scale,tmp,0.0,tmp,Ls-1,Ls-1); } - void PartialFractionFermion5D::M_internal(const LatticeFermion &psi, LatticeFermion &chi,int dag) + template + void PartialFractionFermion5D::M_internal(const FermionField &psi, FermionField &chi,int dag) { - LatticeFermion D(psi._grid); + FermionField D(psi._grid); + int Ls = this->Ls; int sign = dag ? (-1) : 1; // For partial frac Hw case (b5=c5=1) chroma quirkily computes @@ -186,7 +196,7 @@ namespace Grid { // ( 0 -sqrt(p_i)*amax | 2 R gamma_5 + p0/amax 2H // - DW(psi,D,DaggerNo); + this->DW(psi,D,DaggerNo); int nblock=(Ls-1)/2; for(int b=0;b + RealD PartialFractionFermion5D::M (const FermionField &in, FermionField &out) { M_internal(in,out,DaggerNo); return norm2(out); } - RealD PartialFractionFermion5D::Mdag (const LatticeFermion &in, LatticeFermion &out) + template + RealD PartialFractionFermion5D::Mdag (const FermionField &in, FermionField &out) { M_internal(in,out,DaggerYes); return norm2(out); } - void PartialFractionFermion5D::Meooe (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::Meooe (const FermionField &in, FermionField &out) { Meooe_internal(in,out,DaggerNo); } - void PartialFractionFermion5D::MeooeDag (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::MeooeDag (const FermionField &in, FermionField &out) { Meooe_internal(in,out,DaggerYes); } - void PartialFractionFermion5D::Mooee (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::Mooee (const FermionField &in, FermionField &out) { Mooee_internal(in,out,DaggerNo); } - void PartialFractionFermion5D::MooeeDag (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::MooeeDag (const FermionField &in, FermionField &out) { Mooee_internal(in,out,DaggerYes); } - void PartialFractionFermion5D::MooeeInv (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::MooeeInv (const FermionField &in, FermionField &out) { MooeeInv_internal(in,out,DaggerNo); } - void PartialFractionFermion5D::MooeeInvDag (const LatticeFermion &in, LatticeFermion &out) + template + void PartialFractionFermion5D::MooeeInvDag (const FermionField &in, FermionField &out) { MooeeInv_internal(in,out,DaggerYes); } // force terms; five routines; default to Dhop on diagonal - void PartialFractionFermion5D::MDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void PartialFractionFermion5D::MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int nblock=(Ls-1)/2; for(int b=0;bDhopDeriv(mat,D,V,DaggerNo); }; - void PartialFractionFermion5D::MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void PartialFractionFermion5D::MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int nblock=(Ls-1)/2; for(int b=0;bDhopDerivOE(mat,D,V,DaggerNo); }; - void PartialFractionFermion5D::MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) + template + void PartialFractionFermion5D::MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) { - LatticeFermion D(V._grid); + int Ls = this->Ls; + + FermionField D(V._grid); int nblock=(Ls-1)/2; for(int b=0;bDhopDerivEO(mat,D,V,DaggerNo); }; - void PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){ + template + void PartialFractionFermion5D::SetCoefficientsTanh(Approx::zolotarev_data *zdata,RealD scale){ SetCoefficientsZolotarev(1.0/scale,zdata); } - void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){ + template + void PartialFractionFermion5D::SetCoefficientsZolotarev(RealD zolo_hi,Approx::zolotarev_data *zdata){ // check on degree matching // std::cout<db << " -db"<dn << " -dn"<dd << " -dd"<Ls; + assert(Ls == (2*zdata->da -1) ); // Part frac // RealD R; R=(1+mass)/(1-mass); - dw_diag = (4.0-M5); + dw_diag = (4.0-this->M5); // std::vector p; // std::vector q; @@ -336,18 +367,21 @@ namespace Grid { } // Constructors - PartialFractionFermion5D::PartialFractionFermion5D(LatticeGaugeField &_Umu, + template + PartialFractionFermion5D::PartialFractionFermion5D(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, GridRedBlackCartesian &FourDimRedBlackGrid, RealD _mass,RealD M5) : - WilsonFermion5D(_Umu, + WilsonFermion5D(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, FourDimRedBlackGrid,M5), mass(_mass) { + int Ls = this->Ls; + assert((Ls&0x1)==1); // Odd Ls required int nrational=Ls-1; @@ -366,6 +400,8 @@ namespace Grid { } + FermOpTemplateInstantiate(PartialFractionFermion5D); + } } diff --git a/lib/qcd/action/fermion/PartialFractionFermion5D.h b/lib/qcd/action/fermion/PartialFractionFermion5D.h index 3ba2481a..acbc30e1 100644 --- a/lib/qcd/action/fermion/PartialFractionFermion5D.h +++ b/lib/qcd/action/fermion/PartialFractionFermion5D.h @@ -5,46 +5,49 @@ namespace Grid { namespace QCD { - class PartialFractionFermion5D : public WilsonFermion5D + template + class PartialFractionFermion5D : public WilsonFermion5D { public: +#include + public: const int part_frac_chroma_convention=1; - void Meooe_internal(const LatticeFermion &in, LatticeFermion &out,int dag); - void Mooee_internal(const LatticeFermion &in, LatticeFermion &out,int dag); - void MooeeInv_internal(const LatticeFermion &in, LatticeFermion &out,int dag); - void M_internal(const LatticeFermion &in, LatticeFermion &out,int dag); + void Meooe_internal(const FermionField &in, FermionField &out,int dag); + void Mooee_internal(const FermionField &in, FermionField &out,int dag); + void MooeeInv_internal(const FermionField &in, FermionField &out,int dag); + void M_internal(const FermionField &in, FermionField &out,int dag); // override multiply - virtual RealD M (const LatticeFermion &in, LatticeFermion &out); - virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); + virtual RealD M (const FermionField &in, FermionField &out); + virtual RealD Mdag (const FermionField &in, FermionField &out); // half checkerboard operaions - virtual void Meooe (const LatticeFermion &in, LatticeFermion &out); - virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out); - virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); - 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 Meooe (const FermionField &in, FermionField &out); + virtual void MeooeDag (const FermionField &in, FermionField &out); + virtual void Mooee (const FermionField &in, FermionField &out); + virtual void MooeeDag (const FermionField &in, FermionField &out); + virtual void MooeeInv (const FermionField &in, FermionField &out); + virtual void MooeeInvDag (const FermionField &in, FermionField &out); // force terms; five routines; default to Dhop on diagonal - virtual void MDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void MoeDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void MeoDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); virtual void Instantiatable(void) =0; // ensure no make-eee // Efficient support for multigrid coarsening - virtual void Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp); + virtual void Mdir (const FermionField &in, FermionField &out,int dir,int disp); // Constructors - PartialFractionFermion5D(LatticeGaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _mass,RealD M5); + PartialFractionFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _mass,RealD M5); protected: diff --git a/lib/qcd/action/fermion/ScaledShamirFermion.h b/lib/qcd/action/fermion/ScaledShamirFermion.h index 59fb16a8..de7ded0c 100644 --- a/lib/qcd/action/fermion/ScaledShamirFermion.h +++ b/lib/qcd/action/fermion/ScaledShamirFermion.h @@ -7,12 +7,15 @@ namespace Grid { namespace QCD { - class ScaledShamirFermion : public MobiusFermion + template + class ScaledShamirFermion : public MobiusFermion { public: +#include + public: // Constructors - ScaledShamirFermion(LatticeGaugeField &_Umu, + ScaledShamirFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -21,7 +24,7 @@ namespace Grid { RealD scale) : // b+c=scale, b-c = 1 <=> 2b = scale+1; 2c = scale-1 - MobiusFermion(_Umu, + MobiusFermion(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, diff --git a/lib/qcd/action/fermion/ShamirZolotarevFermion.h b/lib/qcd/action/fermion/ShamirZolotarevFermion.h index 6a7df439..339a37f1 100644 --- a/lib/qcd/action/fermion/ShamirZolotarevFermion.h +++ b/lib/qcd/action/fermion/ShamirZolotarevFermion.h @@ -7,14 +7,17 @@ namespace Grid { namespace QCD { - class ShamirZolotarevFermion : public MobiusZolotarevFermion + template + class ShamirZolotarevFermion : public MobiusZolotarevFermion { public: +#include + public: // Constructors - ShamirZolotarevFermion(LatticeGaugeField &_Umu, + ShamirZolotarevFermion(GaugeField &_Umu, GridCartesian &FiveDimGrid, GridRedBlackCartesian &FiveDimRedBlackGrid, GridCartesian &FourDimGrid, @@ -23,7 +26,7 @@ namespace Grid { RealD lo, RealD hi) : // b+c = 1; b-c = 1 => b=1, c=0 - MobiusZolotarevFermion(_Umu, + MobiusZolotarevFermion(_Umu, FiveDimGrid, FiveDimRedBlackGrid, FourDimGrid, diff --git a/lib/qcd/action/fermion/WilsonCompressor.h b/lib/qcd/action/fermion/WilsonCompressor.h index d9fe977a..c7a20c59 100644 --- a/lib/qcd/action/fermion/WilsonCompressor.h +++ b/lib/qcd/action/fermion/WilsonCompressor.h @@ -4,6 +4,7 @@ namespace Grid { namespace QCD { + template class WilsonCompressor { public: int mu; @@ -18,9 +19,13 @@ namespace QCD { mu=p; }; - vHalfSpinColourVector operator () (const vSpinColourVector &in) + virtual SiteHalfSpinor operator () (const SiteSpinor &in) { + return spinproject(in); + } + + SiteHalfSpinor spinproject(const SiteSpinor &in) { - vHalfSpinColourVector ret; + SiteHalfSpinor ret; int mudag=mu; if (dag) { mudag=(mu+Nd)%(2*Nd); @@ -57,5 +62,33 @@ namespace QCD { return ret; } }; + + + template + class GparityWilsonCompressor : public WilsonCompressor{ + public: + + GparityWilsonCompressor(int _dag) : WilsonCompressor (_dag){}; + + SiteHalfSpinor operator () (const SiteSpinor &in) + { + SiteHalfSpinor tmp = spinproject(in); + if( 0 ) tmp = flavourflip(tmp); + return tmp; + } + + SiteHalfSpinor flavourflip(const SiteHalfSpinor &in) { + + SiteHalfSpinor ret; + for(int f=0;f WilsonFermion::directions ({0,1,2,3, 0, 1, 2, 3}); -const std::vector WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1}); + const std::vector WilsonFermionStatic::directions ({0,1,2,3, 0, 1, 2, 3}); + const std::vector WilsonFermionStatic::displacements({1,1,1,1,-1,-1,-1,-1}); + int WilsonFermionStatic::HandOptDslash; -int WilsonFermion::HandOptDslash; + ///////////////////////////////// + // Constructor and gauge import + ///////////////////////////////// -WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu, - GridCartesian &Fgrid, - GridRedBlackCartesian &Hgrid, - RealD _mass) : - _grid(&Fgrid), - _cbgrid(&Hgrid), - Stencil (&Fgrid,npoint,Even,directions,displacements), - StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even - StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd - mass(_mass), - Umu(&Fgrid), - UmuEven(&Hgrid), - UmuOdd (&Hgrid) -{ - // Allocate the required comms buffer - comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO - ImportGauge(_Umu); -} -void WilsonFermion::ImportGauge(const LatticeGaugeField &_Umu) -{ - DoubleStore(Umu,_Umu); - pickCheckerboard(Even,UmuEven,Umu); - pickCheckerboard(Odd ,UmuOdd,Umu); -} - -void WilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) -{ - conformable(Uds._grid,GaugeGrid()); - conformable(Umu._grid,GaugeGrid()); - LatticeColourMatrix U(GaugeGrid()); - for(int mu=0;mu(Umu,mu); - PokeIndex(Uds,U,mu); - U = adj(Cshift(U,mu,-1)); - PokeIndex(Uds,U,mu+4); + template + WilsonFermion::WilsonFermion(GaugeField &_Umu, + GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, + RealD _mass) : + _grid(&Fgrid), + _cbgrid(&Hgrid), + Stencil (&Fgrid,npoint,Even,directions,displacements), + StencilEven(&Hgrid,npoint,Even,directions,displacements), // source is Even + StencilOdd (&Hgrid,npoint,Odd ,directions,displacements), // source is Odd + mass(_mass), + Umu(&Fgrid), + UmuEven(&Hgrid), + UmuOdd (&Hgrid) + { + // Allocate the required comms buffer + comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO + ImportGauge(_Umu); } -} -RealD WilsonFermion::M(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard=in.checkerboard; - Dhop(in,out,DaggerNo); - return axpy_norm(out,4+mass,in,out); -} -RealD WilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard=in.checkerboard; - Dhop(in,out,DaggerYes); - return axpy_norm(out,4+mass,in,out); -} - -void WilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out) -{ - if ( in.checkerboard == Odd ) { - DhopEO(in,out,DaggerNo); - } else { - DhopOE(in,out,DaggerNo); + template + void WilsonFermion::ImportGauge(const GaugeField &_Umu) + { + Impl::DoubleStore(GaugeGrid(),Umu,_Umu); + pickCheckerboard(Even,UmuEven,Umu); + pickCheckerboard(Odd ,UmuOdd,Umu); } -} -void WilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out) -{ - if ( in.checkerboard == Odd ) { - DhopEO(in,out,DaggerYes); - } else { - DhopOE(in,out,DaggerYes); - } -} -void WilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - out = (4.0+mass)*in; - return ; -} -void WilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - Mooee(in,out); -} -void WilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - out = (1.0/(4.0+mass))*in; - return ; -} -void WilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out) -{ - out.checkerboard = in.checkerboard; - MooeeInv(in,out); -} -void WilsonFermion::Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp) -{ - DhopDir(in,out,dir,disp); -} -void WilsonFermion::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp){ - - WilsonCompressor compressor(DaggerNo); - - Stencil.HaloExchange(in,comm_buf,compressor); - assert( (disp==1)||(disp==-1) ); - - int skip = (disp==1) ? 0 : 1; - - int dirdisp = dir+skip*4; - -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,dirdisp); - } - -}; -void WilsonFermion::DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag) -{ - - WilsonCompressor compressor(dag); - - Stencil.HaloExchange(in,comm_buf,compressor); - -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,gamma); - } - -}; - -void WilsonFermion::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, - const LatticeFermion &in, LatticeFermion &out,int dag) -{ - assert((dag==DaggerNo) ||(dag==DaggerYes)); - WilsonCompressor compressor(dag); - st.HaloExchange(in,comm_buf,compressor); - - if ( dag == DaggerYes ) { - if( HandOptDslash ) { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptHandDhopSiteDag(st,U,comm_buf,sss,sss,in,out); - } - } else { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptDhopSiteDag(st,U,comm_buf,sss,sss,in,out); - } - } - } else { - if( HandOptDslash ) { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptHandDhopSite(st,U,comm_buf,sss,sss,in,out); - } - } else { -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptDhopSite(st,U,comm_buf,sss,sss,in,out); - } - } - } -} -void WilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_cbgrid); // verifies half grid - conformable(in._grid,out._grid); // drops the cb check - - assert(in.checkerboard==Even); - out.checkerboard = Odd; - - DhopInternal(StencilEven,UmuOdd,in,out,dag); -} -void WilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_cbgrid); // verifies half grid - conformable(in._grid,out._grid); // drops the cb check - - assert(in.checkerboard==Odd); - out.checkerboard = Even; - - DhopInternal(StencilOdd,UmuEven,in,out,dag); -} -void WilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) -{ - conformable(in._grid,_grid); // verifies full grid - conformable(in._grid,out._grid); - - out.checkerboard = in.checkerboard; - - DhopInternal(Stencil,Umu,in,out,dag); -} - -void WilsonFermion::DerivInternal(CartesianStencil & st,LatticeDoubledGaugeField & U, - LatticeGaugeField &mat,const LatticeFermion &A,const LatticeFermion &B,int dag) -{ - assert((dag==DaggerNo) ||(dag==DaggerYes)); - - WilsonCompressor compressor(dag); - - LatticeColourMatrix tmp(B._grid); - LatticeFermion Btilde(B._grid); - - st.HaloExchange(B,comm_buf,compressor); - - for(int mu=0;mu(1-g) if dag - //////////////////////////////////////////////////////////////////////// - int gamma = mu; - if ( dag ) gamma+= Nd; - - //////////////////////// - // Call the single hop - //////////////////////// -PARALLEL_FOR_LOOP - for(int sss=0;sssoSites();sss++){ - DiracOptDhopDir(st,U,comm_buf,sss,sss,B,Btilde,mu,gamma); - } - - ////////////////////////////////////////////////// - // spin trace outer product - ////////////////////////////////////////////////// - tmp = TraceIndex(outerProduct(Btilde,A)); - PokeIndex(mat,tmp,mu); - + template + RealD WilsonFermion::M(const FermionField &in, FermionField &out) + { + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerNo); + return axpy_norm(out,4+mass,in,out); } -} -void WilsonFermion::DhopDeriv(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) -{ - conformable(U._grid,_grid); - conformable(U._grid,V._grid); - conformable(U._grid,mat._grid); - mat.checkerboard = U.checkerboard; + template + RealD WilsonFermion::Mdag(const FermionField &in, FermionField &out) + { + out.checkerboard=in.checkerboard; + Dhop(in,out,DaggerYes); + return axpy_norm(out,4+mass,in,out); + } - DerivInternal(Stencil,Umu,mat,U,V,dag); -} -void WilsonFermion::DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) -{ - conformable(U._grid,_cbgrid); - conformable(U._grid,V._grid); - conformable(U._grid,mat._grid); + template + void WilsonFermion::Meooe(const FermionField &in, FermionField &out) + { + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerNo); + } else { + DhopOE(in,out,DaggerNo); + } + } + template + void WilsonFermion::MeooeDag(const FermionField &in, FermionField &out) + { + if ( in.checkerboard == Odd ) { + DhopEO(in,out,DaggerYes); + } else { + DhopOE(in,out,DaggerYes); + } + } - assert(V.checkerboard==Even); - assert(U.checkerboard==Odd); - mat.checkerboard = Odd; + template + void WilsonFermion::Mooee(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + out = (4.0+mass)*in; + } + + template + void WilsonFermion::MooeeDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + Mooee(in,out); + } + + template + void WilsonFermion::MooeeInv(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + out = (1.0/(4.0+mass))*in; + } + + template + void WilsonFermion::MooeeInvDag(const FermionField &in, FermionField &out) { + out.checkerboard = in.checkerboard; + MooeeInv(in,out); + } + + /////////////////////////////////// + // Internal + /////////////////////////////////// - DerivInternal(StencilEven,UmuOdd,mat,U,V,dag); -} -void WilsonFermion::DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag) -{ - conformable(U._grid,_cbgrid); - conformable(U._grid,V._grid); - conformable(U._grid,mat._grid); + template + void WilsonFermion::DerivInternal(CartesianStencil & st, + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B,int dag) { + + assert((dag==DaggerNo) ||(dag==DaggerYes)); + + Compressor compressor(dag); + + GaugeLinkField tmp(B._grid); + FermionField Btilde(B._grid); + + st.HaloExchange(B,comm_buf,compressor); + + for(int mu=0;mu(1-g) if dag + //////////////////////////////////////////////////////////////////////// + int gamma = mu; + if ( dag ) gamma+= Nd; + + //////////////////////// + // Call the single hop + //////////////////////// +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptDhopDir(st,U,comm_buf,sss,sss,B,Btilde,mu,gamma); + } + + ////////////////////////////////////////////////// + // spin trace outer product + ////////////////////////////////////////////////// + tmp = TraceIndex(outerProduct(Btilde,A)); + PokeIndex(mat,tmp,mu); + + } + } + + template + void WilsonFermion::DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) + { + conformable(U._grid,_grid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + mat.checkerboard = U.checkerboard; + + DerivInternal(Stencil,Umu,mat,U,V,dag); + } + + template + void WilsonFermion::DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) + { + conformable(U._grid,_cbgrid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + assert(V.checkerboard==Even); + assert(U.checkerboard==Odd); + mat.checkerboard = Odd; + + DerivInternal(StencilEven,UmuOdd,mat,U,V,dag); + } + + template + void WilsonFermion::DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag) + { + conformable(U._grid,_cbgrid); + conformable(U._grid,V._grid); + conformable(U._grid,mat._grid); + + assert(V.checkerboard==Odd); + assert(U.checkerboard==Even); + mat.checkerboard = Even; + + DerivInternal(StencilOdd,UmuEven,mat,U,V,dag); + } + - assert(V.checkerboard==Odd); - assert(U.checkerboard==Even); - mat.checkerboard = Even; + template + void WilsonFermion::Dhop(const FermionField &in, FermionField &out,int dag) { + conformable(in._grid,_grid); // verifies full grid + conformable(in._grid,out._grid); + + out.checkerboard = in.checkerboard; + + DhopInternal(Stencil,Umu,in,out,dag); + } + + template + void WilsonFermion::DhopOE(const FermionField &in, FermionField &out,int dag) { + conformable(in._grid,_cbgrid); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Even); + out.checkerboard = Odd; + + DhopInternal(StencilEven,UmuOdd,in,out,dag); + } + + template + void WilsonFermion::DhopEO(const FermionField &in, FermionField &out,int dag) { + conformable(in._grid,_cbgrid); // verifies half grid + conformable(in._grid,out._grid); // drops the cb check + + assert(in.checkerboard==Odd); + out.checkerboard = Even; + + DhopInternal(StencilOdd,UmuEven,in,out,dag); + } + + template + void WilsonFermion::Mdir (const FermionField &in, FermionField &out,int dir,int disp) { + DhopDir(in,out,dir,disp); + } + + template + void WilsonFermion::DhopDir(const FermionField &in, FermionField &out,int dir,int disp){ + + int skip = (disp==1) ? 0 : 1; + int dirdisp = dir+skip*4; + + DhopDirDisp(in,out,dirdisp,dirdisp,DaggerNo); + + }; + + template + void WilsonFermion::DhopDirDisp(const FermionField &in, FermionField &out,int dirdisp,int gamma,int dag) { + + Compressor compressor(dag); + + Stencil.HaloExchange(in,comm_buf,compressor); + +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptDhopDir(Stencil,Umu,comm_buf,sss,sss,in,out,dirdisp,gamma); + } + + }; - DerivInternal(StencilOdd,UmuEven,mat,U,V,dag); -} + + template + void WilsonFermion::DhopInternal(CartesianStencil & st,DoubledGaugeField & U, + const FermionField &in, FermionField &out,int dag) { + + assert((dag==DaggerNo) ||(dag==DaggerYes)); + + Compressor compressor(dag); + st.HaloExchange(in,comm_buf,compressor); + + if ( dag == DaggerYes ) { + if( HandOptDslash ) { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptHandDhopSiteDag(st,U,comm_buf,sss,sss,in,out); + } + } else { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptDhopSiteDag(st,U,comm_buf,sss,sss,in,out); + } + } + } else { + if( HandOptDslash ) { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptHandDhopSite(st,U,comm_buf,sss,sss,in,out); + } + } else { +PARALLEL_FOR_LOOP + for(int sss=0;sssoSites();sss++){ + Kernels::DiracOptDhopSite(st,U,comm_buf,sss,sss,in,out); + } + } + } + }; + + FermOpTemplateInstantiate(WilsonFermion); }} diff --git a/lib/qcd/action/fermion/WilsonFermion.h b/lib/qcd/action/fermion/WilsonFermion.h index 79d3e626..1dc71c00 100644 --- a/lib/qcd/action/fermion/WilsonFermion.h +++ b/lib/qcd/action/fermion/WilsonFermion.h @@ -5,8 +5,20 @@ namespace Grid { namespace QCD { - class WilsonFermion : public FermionOperator + class WilsonFermionStatic { + public: + static int HandOptDslash; // these are a temporary hack + static int MortonOrder; + static const std::vector directions ; + static const std::vector displacements; + static const int npoint=8; + }; + + template + class WilsonFermion : public FermionOperator, public WilsonFermionStatic { +#include + public: /////////////////////////////////////////////////////////////// @@ -17,133 +29,73 @@ namespace Grid { GridBase *FermionGrid(void) { return _grid;} GridBase *FermionRedBlackGrid(void) { return _cbgrid;} - // override multiply - virtual RealD M (const LatticeFermion &in, LatticeFermion &out); - virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out); + ////////////////////////////////////////////////////////////////// + // override multiply; cut number routines if pass dagger argument + // and also make interface more uniformly consistent + ////////////////////////////////////////////////////////////////// + RealD M(const FermionField &in, FermionField &out); + RealD Mdag(const FermionField &in, FermionField &out); - // half checkerboard operaions - void Meooe (const LatticeFermion &in, LatticeFermion &out); - void MeooeDag (const LatticeFermion &in, LatticeFermion &out); - - virtual void Mooee (const LatticeFermion &in, LatticeFermion &out); // remain virtual so we - virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out); // can derive Clover - virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out); // from Wilson base - virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); + ///////////////////////////////////////////////////////// + // half checkerboard operations + // could remain virtual so we can derive Clover from Wilson base + ///////////////////////////////////////////////////////// + void Meooe(const FermionField &in, FermionField &out) ; + void MeooeDag(const FermionField &in, FermionField &out) ; + void Mooee(const FermionField &in, FermionField &out) ; + void MooeeDag(const FermionField &in, FermionField &out) ; + void MooeeInv(const FermionField &in, FermionField &out) ; + void MooeeInvDag(const FermionField &in, FermionField &out) ; //////////////////////// - // - // Force term: d/dtau S = 0 - // - // It is simplest to consider the two flavour force term - // - // S[U,phi] = phidag (MdagM)^-1 phi - // - // But simplify even this to - // - // S[U,phi] = phidag MdagM phi - // - // (other options exist depending on nature of action fragment.) - // - // Require momentum be traceless anti-hermitian to move within group manifold [ P = i P^a T^a ] - // - // Define the HMC hamiltonian - // - // H = 1/2 Tr P^2 + S(U,phi) - // - // . - // U = P U (lorentz & color indices multiplied) - // - // Hence - // - // .c c c c - // U = U P = - U P (c == dagger) - // - // So, taking some liberty with implicit indices - // . . .c c - // dH/dt = 0 = Tr P P +Tr[ U dS/dU + U dS/dU ] - // - // . c c - // = Tr P P + i Tr[ P U dS/dU - U P dS/dU ] - // - // . c c - // = Tr P (P + i ( U dS/dU - P dS/dU U ] - // - // . c c - // => P = -i [ U dS/dU - dS/dU U ] generates HMC EoM - // - // Simple case work this out using S = phi^dag MdagM phi for wilson: - // c c - // dSdt = dU_xdt dSdUx + dUxdt dSdUx - // - // = Tr i P U_x [ (\phi^\dag)_x (1+g) (M \phi)_x+\mu +(\phi^\dag M^\dag)_x (1-g) \phi_{x+\mu} ] - // c - // - i U_x P [ (\phi^\dag)_x+mu (1-g) (M \phi)_x +(\phi^\dag M^\dag)_(x+\mu) (1+g) \phi_{x} ] - // - // = i [(\phi^\dag)_x ]_j P_jk [U_x(1+g) (M \phi)_x+\mu]_k (1) - // + i [(\phi^\dagM^\dag)_x]_j P_jk [U_x(1-g) (\phi)_x+\mu]_k (2) - // - i [(\phi^\dag)_x+mu (1-g) U^dag_x]_j P_jk [(M \phi)_xk (3) - // - i [(\phi^\dagM^\dag)_x+mu (1+g) U^dag_x]_j P_jk [ \phi]_xk (4) - // - // Observe that (1)* = (4) - // (2)* = (3) - // - // Write as . - // P_{kj} = - i ( [U_x(1+g) (M \phi)_x+\mu] (x) [(\phi^\dag)_x] + [U_x(1-g) (\phi)_x+\mu] (x) [(\phi^\dagM^\dag)_x] - h.c ) - // - // where (x) denotes outer product in colour and spins are traced. - // - // Need only evaluate (1) and (2) [Chroma] or (2) and (4) [IroIro] and take the - // traceless anti hermitian part (of term in brackets w/o the "i") - // - // Generalisation to S=phi^dag (MdagM)^{-1} phi is simple: - // - // For more complicated DWF etc... apply product rule in differentiation - // + // Derivative interface //////////////////////// - void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - - // Extra support internal - void DerivInternal(CartesianStencil & st, - LatticeDoubledGaugeField & U, - LatticeGaugeField &mat, - const LatticeFermion &A, - const LatticeFermion &B, - int dag); + // Interface calls an internal routine + void DhopDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + /////////////////////////////////////////////////////////////// // non-hermitian hopping term; half cb or both - void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag); + /////////////////////////////////////////////////////////////// + void Dhop(const FermionField &in, FermionField &out,int dag) ; + void DhopOE(const FermionField &in, FermionField &out,int dag) ; + void DhopEO(const FermionField &in, FermionField &out,int dag) ; - // Multigrid assistance - void Mdir (const LatticeFermion &in, LatticeFermion &out,int dir,int disp); - void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp); - void DhopDirDisp(const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma,int dag); + /////////////////////////////////////////////////////////////// + // Multigrid assistance; force term uses too + /////////////////////////////////////////////////////////////// + void Mdir (const FermionField &in, FermionField &out,int dir,int disp) ; + void DhopDir(const FermionField &in, FermionField &out,int dir,int disp); + void DhopDirDisp(const FermionField &in, FermionField &out,int dirdisp,int gamma,int dag) ; /////////////////////////////////////////////////////////////// // Extra methods added by derived /////////////////////////////////////////////////////////////// - void DhopInternal(CartesianStencil & st, - LatticeDoubledGaugeField &U, - const LatticeFermion &in, - LatticeFermion &out, - int dag); + void DerivInternal(CartesianStencil & st, + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag); + + void DhopInternal(CartesianStencil & st,DoubledGaugeField & U, + const FermionField &in, FermionField &out,int dag) ; + // Constructor - WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass); + WilsonFermion(GaugeField &_Umu, + GridCartesian &Fgrid, + GridRedBlackCartesian &Hgrid, + RealD _mass) ; - // DoubleStore - virtual void ImportGauge(const LatticeGaugeField &_Umu); - void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); + // DoubleStore impl dependent + void ImportGauge(const GaugeField &_Umu); /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// - static int HandOptDslash; // these are a temporary hack - static int MortonOrder; // protected: public: @@ -153,26 +105,24 @@ namespace Grid { GridBase * _grid; GridBase * _cbgrid; - static const int npoint=8; - static const std::vector directions ; - static const std::vector displacements; - //Defines the stencils for even and odd CartesianStencil Stencil; CartesianStencil StencilEven; CartesianStencil StencilOdd; // Copy of the gauge field , with even and odd subsets - LatticeDoubledGaugeField Umu; - LatticeDoubledGaugeField UmuEven; - LatticeDoubledGaugeField UmuOdd; + DoubledGaugeField Umu; + DoubledGaugeField UmuEven; + DoubledGaugeField UmuOdd; // Comms buffer - std::vector > comm_buf; - + std::vector > comm_buf; }; + typedef WilsonFermion WilsonFermionF; + typedef WilsonFermion WilsonFermionD; + } } #endif diff --git a/lib/qcd/action/fermion/WilsonFermion5D.cc b/lib/qcd/action/fermion/WilsonFermion5D.cc index cbe66bb5..f324db8d 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.cc +++ b/lib/qcd/action/fermion/WilsonFermion5D.cc @@ -3,19 +3,19 @@ namespace Grid { namespace QCD { - // S-direction is INNERMOST and takes no part in the parity. - const std::vector WilsonFermion5D::directions ({1,2,3,4, 1, 2, 3, 4}); - const std::vector WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1}); - - int WilsonFermion5D::HandOptDslash; +// S-direction is INNERMOST and takes no part in the parity. +const std::vector WilsonFermion5DStatic::directions ({1,2,3,4, 1, 2, 3, 4}); +const std::vector WilsonFermion5DStatic::displacements({1,1,1,1,-1,-1,-1,-1}); +int WilsonFermion5DStatic::HandOptDslash; // 5d lattice for DWF. - WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - RealD _M5) : +template +WilsonFermion5D::WilsonFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + RealD _M5) : _FiveDimGrid(&FiveDimGrid), _FiveDimRedBlackGrid(&FiveDimRedBlackGrid), _FourDimGrid(&FourDimGrid), @@ -68,34 +68,23 @@ namespace QCD { ImportGauge(_Umu); } -void WilsonFermion5D::ImportGauge(const LatticeGaugeField &_Umu) +template +void WilsonFermion5D::ImportGauge(const GaugeField &_Umu) { - DoubleStore(Umu,_Umu); + Impl::DoubleStore(GaugeGrid(),Umu,_Umu); pickCheckerboard(Even,UmuEven,Umu); pickCheckerboard(Odd ,UmuOdd,Umu); } -void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) -{ - assert(GaugeGrid()->_ndimension==4); - conformable(Uds._grid,GaugeGrid()); - conformable(Umu._grid,GaugeGrid()); - LatticeColourMatrix U(GaugeGrid()); - for(int mu=0;mu(Umu,mu); - PokeIndex(Uds,U,mu); - U = adj(Cshift(U,mu,-1)); - PokeIndex(Uds,U,mu+4); - } -} -void WilsonFermion5D::DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir5,int disp) +template +void WilsonFermion5D::DhopDir(const FermionField &in, FermionField &out,int dir5,int disp) { int dir = dir5-1; // Maps to the ordering above in "directions" that is passed to stencil // we drop off the innermost fifth dimension // assert( (disp==1)||(disp==-1) ); // assert( (dir>=0)&&(dir<4) ); //must do x,y,z or t; - WilsonCompressor compressor(DaggerNo); - Stencil.HaloExchange(in,comm_buf,compressor); + Compressor compressor(DaggerNo); + Stencil.HaloExchange(in,comm_buf,compressor); int skip = (disp==1) ? 0 : 1; @@ -109,16 +98,17 @@ PARALLEL_FOR_LOOP for(int s=0;s +void WilsonFermion5D::DerivInternal(CartesianStencil & st, + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B, int dag) { assert((dag==DaggerNo) ||(dag==DaggerYes)); @@ -126,13 +116,13 @@ void WilsonFermion5D::DerivInternal(CartesianStencil & st, conformable(st._grid,A._grid); conformable(st._grid,B._grid); - WilsonCompressor compressor(dag); + Compressor compressor(dag); - LatticeColourMatrix tmp(mat._grid); - LatticeFermion Btilde(B._grid); - LatticeFermion Atilde(B._grid); + GaugeLinkField tmp(mat._grid); + FermionField Btilde(B._grid); + FermionField Atilde(B._grid); - st.HaloExchange(B,comm_buf,compressor); + st.HaloExchange(B,comm_buf,compressor); Atilde=A; @@ -158,7 +148,7 @@ PARALLEL_FOR_LOOP assert ( sF< B._grid->oSites()); assert ( sU< U._grid->oSites()); - DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma); + Kernels::DiracOptDhopDir(st,U,comm_buf,sF,sU,B,Btilde,mu,gamma); //////////////////////////// // spin trace outer product @@ -176,10 +166,11 @@ PARALLEL_FOR_LOOP } } -void WilsonFermion5D::DhopDeriv( LatticeGaugeField &mat, - const LatticeFermion &A, - const LatticeFermion &B, - int dag) +template +void WilsonFermion5D::DhopDeriv( GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) { conformable(A._grid,FermionGrid()); conformable(A._grid,B._grid); @@ -190,10 +181,11 @@ void WilsonFermion5D::DhopDeriv( LatticeGaugeField &mat, DerivInternal(Stencil,Umu,mat,A,B,dag); } -void WilsonFermion5D::DhopDerivEO(LatticeGaugeField &mat, - const LatticeFermion &A, - const LatticeFermion &B, - int dag) +template +void WilsonFermion5D::DhopDerivEO(GaugeField &mat, + const FermionField &A, + const FermionField &B, + int dag) { conformable(A._grid,FermionRedBlackGrid()); conformable(GaugeRedBlackGrid(),mat._grid); @@ -206,9 +198,10 @@ void WilsonFermion5D::DhopDerivEO(LatticeGaugeField &mat, DerivInternal(StencilOdd,UmuEven,mat,A,B,dag); } -void WilsonFermion5D::DhopDerivOE(LatticeGaugeField &mat, - const LatticeFermion &A, - const LatticeFermion &B, +template +void WilsonFermion5D::DhopDerivOE(GaugeField &mat, + const FermionField &A, + const FermionField &B, int dag) { conformable(A._grid,FermionRedBlackGrid()); @@ -222,15 +215,16 @@ void WilsonFermion5D::DhopDerivOE(LatticeGaugeField &mat, DerivInternal(StencilEven,UmuOdd,mat,A,B,dag); } -void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo, - LatticeDoubledGaugeField & U, - const LatticeFermion &in, LatticeFermion &out,int dag) +template +void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo, + DoubledGaugeField & U, + const FermionField &in, FermionField &out,int dag) { // assert((dag==DaggerNo) ||(dag==DaggerYes)); - WilsonCompressor compressor(dag); + Compressor compressor(dag); - st.HaloExchange(in,comm_buf,compressor); + st.HaloExchange(in,comm_buf,compressor); // Dhop takes the 4d grid from U, and makes a 5d index for fermion // Not loop ordering and data layout. @@ -238,13 +232,13 @@ void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo, // - per thread reuse in L1 cache for U // - 8 linear access unit stride streams per thread for Fermion for hw prefetchable. if ( dag == DaggerYes ) { - if( HandOptDslash ) { + if( this->HandOptDslash ) { PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ for(int s=0;sHandOptDslash ) { PARALLEL_FOR_LOOP for(int ss=0;ssoSites();ss++){ for(int s=0;s +void WilsonFermion5D::DhopOE(const FermionField &in, FermionField &out,int dag) { conformable(in._grid,FermionRedBlackGrid()); // verifies half grid conformable(in._grid,out._grid); // drops the cb check @@ -295,7 +290,8 @@ void WilsonFermion5D::DhopOE(const LatticeFermion &in, LatticeFermion &out,int d DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag); } -void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag) +template +void WilsonFermion5D::DhopEO(const FermionField &in, FermionField &out,int dag) { conformable(in._grid,FermionRedBlackGrid()); // verifies half grid conformable(in._grid,out._grid); // drops the cb check @@ -305,7 +301,8 @@ void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int d DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag); } -void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) +template +void WilsonFermion5D::Dhop(const FermionField &in, FermionField &out,int dag) { conformable(in._grid,FermionGrid()); // verifies full grid conformable(in._grid,out._grid); @@ -314,12 +311,16 @@ void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag DhopInternal(Stencil,Lebesgue,Umu,in,out,dag); } -void WilsonFermion5D::DW(const LatticeFermion &in, LatticeFermion &out,int dag) +template +void WilsonFermion5D::DW(const FermionField &in, FermionField &out,int dag) { out.checkerboard=in.checkerboard; Dhop(in,out,dag); // -0.5 is included axpy(out,4.0-M5,in,out); } + +FermOpTemplateInstantiate(WilsonFermion5D); + }} diff --git a/lib/qcd/action/fermion/WilsonFermion5D.h b/lib/qcd/action/fermion/WilsonFermion5D.h index d5cff3e6..5159d9e1 100644 --- a/lib/qcd/action/fermion/WilsonFermion5D.h +++ b/lib/qcd/action/fermion/WilsonFermion5D.h @@ -14,21 +14,21 @@ 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 + + class WilsonFermion5DStatic { + public: + // S-direction is INNERMOST and takes no part in the parity. + static int HandOptDslash; // these are a temporary hack + static const std::vector directions; + static const std::vector displacements; + const int npoint = 8; + }; + + template + class WilsonFermion5D : public FermionOperator, public WilsonFermion5DStatic { +#include public: /////////////////////////////////////////////////////////////// // Implement the abstract base @@ -39,69 +39,65 @@ namespace Grid { GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;} // full checkerboard operations; leave unimplemented as abstract for now - virtual RealD M (const LatticeFermion &in, LatticeFermion &out){assert(0); return 0.0;}; - virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out){assert(0); return 0.0;}; + virtual RealD M (const FermionField &in, FermionField &out){assert(0); return 0.0;}; + virtual RealD Mdag (const FermionField &in, FermionField &out){assert(0); return 0.0;}; // half checkerboard operations; leave unimplemented as abstract for now - virtual void Meooe (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void Mooee (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out){assert(0);}; + virtual void Meooe (const FermionField &in, FermionField &out){assert(0);}; + virtual void Mooee (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeInv (const FermionField &in, FermionField &out){assert(0);}; - virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; - virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out){assert(0);}; + virtual void MeooeDag (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeDag (const FermionField &in, FermionField &out){assert(0);}; + virtual void MooeeInvDag (const FermionField &in, FermionField &out){assert(0);}; - // These can be overridden by fancy 5d chiral actions - virtual void DhopDeriv (LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void DhopDerivEO(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); - virtual void DhopDerivOE(LatticeGaugeField &mat,const LatticeFermion &U,const LatticeFermion &V,int dag); + // These can be overridden by fancy 5d chiral action + virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); + virtual void DhopDerivOE(GaugeField &mat,const FermionField &U,const FermionField &V,int dag); // Implement hopping term non-hermitian hopping term; half cb or both // Implement s-diagonal DW - void DW (const LatticeFermion &in, LatticeFermion &out,int dag); - void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag); - void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag); + void DW (const FermionField &in, FermionField &out,int dag); + void Dhop (const FermionField &in, FermionField &out,int dag); + void DhopOE(const FermionField &in, FermionField &out,int dag); + void DhopEO(const FermionField &in, FermionField &out,int dag); // add a DhopComm // -- suboptimal interface will presently trigger multiple comms. - void DhopDir(const LatticeFermion &in, LatticeFermion &out,int dir,int disp); + void DhopDir(const FermionField &in, FermionField &out,int dir,int disp); /////////////////////////////////////////////////////////////// // New methods added /////////////////////////////////////////////////////////////// - void DerivInternal(CartesianStencil & st, - LatticeDoubledGaugeField & U, - LatticeGaugeField &mat, - const LatticeFermion &A, - const LatticeFermion &B, + DoubledGaugeField & U, + GaugeField &mat, + const FermionField &A, + const FermionField &B, int dag); void DhopInternal(CartesianStencil & st, LebesgueOrder &lo, - LatticeDoubledGaugeField &U, - const LatticeFermion &in, - LatticeFermion &out, + DoubledGaugeField &U, + const FermionField &in, + FermionField &out, int dag); // Constructors - WilsonFermion5D(LatticeGaugeField &_Umu, - GridCartesian &FiveDimGrid, - GridRedBlackCartesian &FiveDimRedBlackGrid, - GridCartesian &FourDimGrid, - GridRedBlackCartesian &FourDimRedBlackGrid, - double _M5); + WilsonFermion5D(GaugeField &_Umu, + GridCartesian &FiveDimGrid, + GridRedBlackCartesian &FiveDimRedBlackGrid, + GridCartesian &FourDimGrid, + GridRedBlackCartesian &FourDimRedBlackGrid, + double _M5); // DoubleStore - virtual void ImportGauge(const LatticeGaugeField &_Umu); - void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); + void ImportGauge(const GaugeField &_Umu); /////////////////////////////////////////////////////////////// // Data members require to support the functionality /////////////////////////////////////////////////////////////// - static int HandOptDslash; // these are a temporary hack - protected: // Add these to the support from Wilson @@ -110,10 +106,6 @@ namespace Grid { GridBase *_FiveDimGrid; GridBase *_FiveDimRedBlackGrid; - static const int npoint=8; - static const std::vector directions ; - static const std::vector displacements; - double M5; int Ls; @@ -123,15 +115,15 @@ namespace Grid { CartesianStencil StencilOdd; // Copy of the gauge field , with even and odd subsets - LatticeDoubledGaugeField Umu; - LatticeDoubledGaugeField UmuEven; - LatticeDoubledGaugeField UmuOdd; + DoubledGaugeField Umu; + DoubledGaugeField UmuEven; + DoubledGaugeField UmuOdd; LebesgueOrder Lebesgue; LebesgueOrder LebesgueEvenOdd; // Comms buffer - std::vector > comm_buf; + std::vector > comm_buf; }; } diff --git a/lib/qcd/action/fermion/WilsonKernels.cc b/lib/qcd/action/fermion/WilsonKernels.cc index 2b8b7929..6d558e03 100644 --- a/lib/qcd/action/fermion/WilsonKernels.cc +++ b/lib/qcd/action/fermion/WilsonKernels.cc @@ -2,23 +2,317 @@ namespace Grid { namespace QCD { -void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out) +template +void WilsonKernels::DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out) { - vHalfSpinColourVector tmp; - vHalfSpinColourVector chi; - vSpinColourVector result; - vHalfSpinColourVector Uchi; - int offset,local,perm, ptype; + SiteHalfSpinor tmp; + SiteHalfSpinor chi; + SiteSpinor result; + SiteHalfSpinor Uchi; + int offset,local,perm, ptype; - // Xp - int ss = sF; - offset = st._offsets [Xp][ss]; - local = st._is_local[Xp][ss]; - perm = st._permute[Xp][ss]; + // Xp + int ss = sF; + offset = st._offsets [Xp][ss]; + local = st._is_local[Xp][ss]; + perm = st._permute[Xp][ss]; + + ptype = st._permute_type[Xp]; + if ( local && perm ) { + spProjXp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjXp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Xp); + spReconXp(result,Uchi); + + // Yp + offset = st._offsets [Yp][ss]; + local = st._is_local[Yp][ss]; + perm = st._permute[Yp][ss]; + ptype = st._permute_type[Yp]; + if ( local && perm ) { + spProjYp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjYp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Yp); + accumReconYp(result,Uchi); - ptype = st._permute_type[Xp]; + // Zp + offset = st._offsets [Zp][ss]; + local = st._is_local[Zp][ss]; + perm = st._permute[Zp][ss]; + ptype = st._permute_type[Zp]; + if ( local && perm ) { + spProjZp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjZp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Zp); + accumReconZp(result,Uchi); + + // Tp + offset = st._offsets [Tp][ss]; + local = st._is_local[Tp][ss]; + perm = st._permute[Tp][ss]; + ptype = st._permute_type[Tp]; + if ( local && perm ) { + spProjTp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjTp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Tp); + accumReconTp(result,Uchi); + + // Xm + offset = st._offsets [Xm][ss]; + local = st._is_local[Xm][ss]; + perm = st._permute[Xm][ss]; + ptype = st._permute_type[Xm]; + + if ( local && perm ) { + spProjXm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjXm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Xm); + accumReconXm(result,Uchi); + + // Ym + offset = st._offsets [Ym][ss]; + local = st._is_local[Ym][ss]; + perm = st._permute[Ym][ss]; + ptype = st._permute_type[Ym]; + + if ( local && perm ) { + spProjYm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjYm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Ym); + accumReconYm(result,Uchi); + + // Zm + offset = st._offsets [Zm][ss]; + local = st._is_local[Zm][ss]; + perm = st._permute[Zm][ss]; + ptype = st._permute_type[Zm]; + if ( local && perm ) { + spProjZm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjZm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Zm); + accumReconZm(result,Uchi); + + // Tm + offset = st._offsets [Tm][ss]; + local = st._is_local[Tm][ss]; + perm = st._permute[Tm][ss]; + ptype = st._permute_type[Tm]; + if ( local && perm ) { + spProjTm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjTm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Tm); + accumReconTm(result,Uchi); + + vstream(out._odata[ss],result*(-0.5)); +}; + +template +void WilsonKernels::DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out) +{ + SiteHalfSpinor tmp; + SiteHalfSpinor chi; + SiteSpinor result; + SiteHalfSpinor Uchi; + int offset,local,perm, ptype; + + // Xp + int ss=sF; + offset = st._offsets [Xm][ss]; + local = st._is_local[Xm][ss]; + perm = st._permute[Xm][ss]; + + ptype = st._permute_type[Xm]; + if ( local && perm ) { + spProjXp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjXp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Xm); + spReconXp(result,Uchi); + + // Yp + offset = st._offsets [Ym][ss]; + local = st._is_local[Ym][ss]; + perm = st._permute[Ym][ss]; + ptype = st._permute_type[Ym]; + if ( local && perm ) { + spProjYp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjYp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Ym); + accumReconYp(result,Uchi); + + // Zp + offset = st._offsets [Zm][ss]; + local = st._is_local[Zm][ss]; + perm = st._permute[Zm][ss]; + ptype = st._permute_type[Zm]; + if ( local && perm ) { + spProjZp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjZp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Zm); + accumReconZp(result,Uchi); + + // Tp + offset = st._offsets [Tm][ss]; + local = st._is_local[Tm][ss]; + perm = st._permute[Tm][ss]; + ptype = st._permute_type[Tm]; + if ( local && perm ) { + spProjTp(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjTp(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Tm); + accumReconTp(result,Uchi); + + // Xm + offset = st._offsets [Xp][ss]; + local = st._is_local[Xp][ss]; + perm = st._permute[Xp][ss]; + ptype = st._permute_type[Xp]; + + if ( local && perm ) { + spProjXm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjXm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Xp); + accumReconXm(result,Uchi); + + // Ym + offset = st._offsets [Yp][ss]; + local = st._is_local[Yp][ss]; + perm = st._permute[Yp][ss]; + ptype = st._permute_type[Yp]; + + if ( local && perm ) { + spProjYm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjYm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Yp); + accumReconYm(result,Uchi); + + // Zm + offset = st._offsets [Zp][ss]; + local = st._is_local[Zp][ss]; + perm = st._permute[Zp][ss]; + ptype = st._permute_type[Zp]; + if ( local && perm ) { + spProjZm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjZm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Zp); + accumReconZm(result,Uchi); + + // Tm + offset = st._offsets [Tp][ss]; + local = st._is_local[Tp][ss]; + perm = st._permute[Tp][ss]; + ptype = st._permute_type[Tp]; + if ( local && perm ) { + spProjTm(tmp,in._odata[offset]); + permute(chi,tmp,ptype); + } else if ( local ) { + spProjTm(chi,in._odata[offset]); + } else { + chi=buf[offset]; + } + Impl::multLink(Uchi,U._odata[sU],chi,Tp); + accumReconTm(result,Uchi); + + vstream(out._odata[ss],result*(-0.5)); +} + +template +void WilsonKernels::DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out,int dir,int gamma) +{ + SiteHalfSpinor tmp; + SiteHalfSpinor chi; + SiteSpinor result; + SiteHalfSpinor Uchi; + int offset,local,perm, ptype; + int ss=sF; + + offset = st._offsets [dir][ss]; + local = st._is_local[dir][ss]; + perm = st._permute[dir][ss]; + ptype = st._permute_type[dir]; + + // Xp + if(gamma==Xp){ if ( local && perm ) { spProjXp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -27,14 +321,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Xp),&chi()); + Impl::multLink(Uchi,U._odata[sU],chi,dir); spReconXp(result,Uchi); + } - // Yp - offset = st._offsets [Yp][ss]; - local = st._is_local[Yp][ss]; - perm = st._permute[Yp][ss]; - ptype = st._permute_type[Yp]; + // Yp + if ( gamma==Yp ){ if ( local && perm ) { spProjYp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -43,14 +335,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Yp),&chi()); - accumReconYp(result,Uchi); - - // Zp - offset = st._offsets [Zp][ss]; - local = st._is_local[Zp][ss]; - perm = st._permute[Zp][ss]; - ptype = st._permute_type[Zp]; + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconYp(result,Uchi); + } + + // Zp + if ( gamma ==Zp ){ if ( local && perm ) { spProjZp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -59,14 +349,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Zp),&chi()); - accumReconZp(result,Uchi); - - // Tp - offset = st._offsets [Tp][ss]; - local = st._is_local[Tp][ss]; - perm = st._permute[Tp][ss]; - ptype = st._permute_type[Tp]; + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconZp(result,Uchi); + } + + // Tp + if ( gamma ==Tp ){ if ( local && perm ) { spProjTp(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -75,15 +363,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Tp),&chi()); - accumReconTp(result,Uchi); - - // Xm - offset = st._offsets [Xm][ss]; - local = st._is_local[Xm][ss]; - perm = st._permute[Xm][ss]; - ptype = st._permute_type[Xm]; + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconTp(result,Uchi); + } + // Xm + if ( gamma==Xm ){ if ( local && perm ) { spProjXm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -92,15 +377,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Xm),&chi()); - accumReconXm(result,Uchi); - - // Ym - offset = st._offsets [Ym][ss]; - local = st._is_local[Ym][ss]; - perm = st._permute[Ym][ss]; - ptype = st._permute_type[Ym]; + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconXm(result,Uchi); + } + // Ym + if ( gamma == Ym ){ if ( local && perm ) { spProjYm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -109,14 +391,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Ym),&chi()); - accumReconYm(result,Uchi); + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconYm(result,Uchi); + } - // Zm - offset = st._offsets [Zm][ss]; - local = st._is_local[Zm][ss]; - perm = st._permute[Zm][ss]; - ptype = st._permute_type[Zm]; + // Zm + if ( gamma == Zm ){ if ( local && perm ) { spProjZm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -125,14 +405,12 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Zm),&chi()); - accumReconZm(result,Uchi); - - // Tm - offset = st._offsets [Tm][ss]; - local = st._is_local[Tm][ss]; - perm = st._permute[Tm][ss]; - ptype = st._permute_type[Tm]; + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconZm(result,Uchi); + } + + // Tm + if ( gamma==Tm ) { if ( local && perm ) { spProjTm(tmp,in._odata[offset]); permute(chi,tmp,ptype); @@ -141,287 +419,13 @@ void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, } else { chi=buf[offset]; } - mult(&Uchi(),&U._odata[sU](Tm),&chi()); - accumReconTm(result,Uchi); + Impl::multLink(Uchi,U._odata[sU],chi,dir); + spReconTm(result,Uchi); + } - vstream(out._odata[ss],result*(-0.5)); + vstream(out._odata[ss],result*(-0.5)); } -void DiracOptDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out) -{ - vHalfSpinColourVector tmp; - vHalfSpinColourVector chi; - vSpinColourVector result; - vHalfSpinColourVector Uchi; - int offset,local,perm, ptype; - - // Xp - int ss=sF; - offset = st._offsets [Xm][ss]; - local = st._is_local[Xm][ss]; - perm = st._permute[Xm][ss]; - - ptype = st._permute_type[Xm]; - if ( local && perm ) { - spProjXp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjXp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Xm),&chi()); - spReconXp(result,Uchi); - - // Yp - offset = st._offsets [Ym][ss]; - local = st._is_local[Ym][ss]; - perm = st._permute[Ym][ss]; - ptype = st._permute_type[Ym]; - if ( local && perm ) { - spProjYp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjYp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Ym),&chi()); - accumReconYp(result,Uchi); - - // Zp - offset = st._offsets [Zm][ss]; - local = st._is_local[Zm][ss]; - perm = st._permute[Zm][ss]; - ptype = st._permute_type[Zm]; - if ( local && perm ) { - spProjZp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjZp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Zm),&chi()); - accumReconZp(result,Uchi); - - // Tp - offset = st._offsets [Tm][ss]; - local = st._is_local[Tm][ss]; - perm = st._permute[Tm][ss]; - ptype = st._permute_type[Tm]; - if ( local && perm ) { - spProjTp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjTp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Tm),&chi()); - accumReconTp(result,Uchi); - - // Xm - offset = st._offsets [Xp][ss]; - local = st._is_local[Xp][ss]; - perm = st._permute[Xp][ss]; - ptype = st._permute_type[Xp]; - - if ( local && perm ) - { - spProjXm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjXm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Xp),&chi()); - accumReconXm(result,Uchi); - - // Ym - offset = st._offsets [Yp][ss]; - local = st._is_local[Yp][ss]; - perm = st._permute[Yp][ss]; - ptype = st._permute_type[Yp]; - - if ( local && perm ) { - spProjYm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjYm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Yp),&chi()); - accumReconYm(result,Uchi); - - // Zm - offset = st._offsets [Zp][ss]; - local = st._is_local[Zp][ss]; - perm = st._permute[Zp][ss]; - ptype = st._permute_type[Zp]; - if ( local && perm ) { - spProjZm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjZm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Zp),&chi()); - accumReconZm(result,Uchi); - - // Tm - offset = st._offsets [Tp][ss]; - local = st._is_local[Tp][ss]; - perm = st._permute[Tp][ss]; - ptype = st._permute_type[Tp]; - if ( local && perm ) { - spProjTm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjTm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](Tp),&chi()); - accumReconTm(result,Uchi); - - vstream(out._odata[ss],result*(-0.5)); -} - -void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dir,int gamma) -{ - vHalfSpinColourVector tmp; - vHalfSpinColourVector chi; - vSpinColourVector result; - vHalfSpinColourVector Uchi; - int offset,local,perm, ptype; - int ss=sF; - - offset = st._offsets [dir][ss]; - local = st._is_local[dir][ss]; - perm = st._permute[dir][ss]; - ptype = st._permute_type[dir]; - - // Xp - if(gamma==Xp){ - if ( local && perm ) { - spProjXp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjXp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconXp(result,Uchi); - } - - // Yp - if ( gamma==Yp ){ - if ( local && perm ) { - spProjYp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjYp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconYp(result,Uchi); - } - - // Zp - if ( gamma ==Zp ){ - if ( local && perm ) { - spProjZp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjZp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconZp(result,Uchi); - } - - // Tp - if ( gamma ==Tp ){ - if ( local && perm ) { - spProjTp(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjTp(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconTp(result,Uchi); - } - - // Xm - if ( gamma==Xm ){ - if ( local && perm ) { - spProjXm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjXm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconXm(result,Uchi); - } - - // Ym - if ( gamma == Ym ){ - if ( local && perm ) { - spProjYm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjYm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconYm(result,Uchi); - } - - // Zm - if ( gamma == Zm ){ - if ( local && perm ) { - spProjZm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjZm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconZm(result,Uchi); - } - - // Tm - if ( gamma==Tm ) { - if ( local && perm ) { - spProjTm(tmp,in._odata[offset]); - permute(chi,tmp,ptype); - } else if ( local ) { - spProjTm(chi,in._odata[offset]); - } else { - chi=buf[offset]; - } - mult(&Uchi(),&U._odata[sU](dir),&chi()); - spReconTm(result,Uchi); - } - - vstream(out._odata[ss],result*(-0.5)); -} + FermOpTemplateInstantiate(WilsonKernels); }} diff --git a/lib/qcd/action/fermion/WilsonKernels.h b/lib/qcd/action/fermion/WilsonKernels.h index cb0e07bb..eb0981fc 100644 --- a/lib/qcd/action/fermion/WilsonKernels.h +++ b/lib/qcd/action/fermion/WilsonKernels.h @@ -6,44 +6,46 @@ namespace Grid { namespace QCD { //////////////////////////////////////////////////////////////////////////////////////////////////////////////// - // Helper classes that implement Wilson stencil for a single site. + // Helper routines that implement Wilson stencil for a single site. //////////////////////////////////////////////////////////////////////////////////////////////////////////////// - - // Generic version works for any Nc and with extra flavour indices - // namespace DiracOpt { - // These ones will need to be package intelligently. WilsonType base class - // for use by DWF etc.. - void DiracOptDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out); - void DiracOptDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out); - void DiracOptDhopDir(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out,int dirdisp,int gamma); + template + class WilsonKernels { + public: + +#include + + public: + static void + DiracOptDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out); - // }; + static void + DiracOptDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in,FermionField &out); - // Hand unrolled for Nc=3, one flavour - // namespace DiracOptHand { - // These ones will need to be package intelligently. WilsonType base class - // for use by DWF etc.. + static void + DiracOptDhopDir(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out,int dirdisp,int gamma); - void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out); - void DiracOptHandDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + static void + DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out){ + DiracOptDhopSite(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 + } - // }; - + static void + DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out){ + DiracOptDhopSiteDag(st,U,buf,sF,sU,in,out); // will template override for Wilson Nc=3 + } - void DiracOptHandDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out); + }; } } diff --git a/lib/qcd/action/fermion/WilsonKernelsHand.cc b/lib/qcd/action/fermion/WilsonKernelsHand.cc index cff700fb..ea9d5457 100644 --- a/lib/qcd/action/fermion/WilsonKernelsHand.cc +++ b/lib/qcd/action/fermion/WilsonKernelsHand.cc @@ -280,48 +280,50 @@ namespace Grid { namespace QCD { -void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, - std::vector > &buf, - int sF,int sU,const LatticeFermion &in, LatticeFermion &out) +#if 0 +template +void WilsonKernels >::DiracOptHandDhopSite(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int sF,int sU,const FermionField &in, FermionField &out) { - REGISTER vComplex result_00; // 12 regs on knc - REGISTER vComplex result_01; - REGISTER vComplex result_02; + REGISTER Simd result_00; // 12 regs on knc + REGISTER Simd result_01; + REGISTER Simd result_02; - REGISTER vComplex result_10; - REGISTER vComplex result_11; - REGISTER vComplex result_12; + REGISTER Simd result_10; + REGISTER Simd result_11; + REGISTER Simd result_12; - REGISTER vComplex result_20; - REGISTER vComplex result_21; - REGISTER vComplex result_22; + REGISTER Simd result_20; + REGISTER Simd result_21; + REGISTER Simd result_22; - REGISTER vComplex result_30; - REGISTER vComplex result_31; - REGISTER vComplex result_32; // 20 left + REGISTER Simd result_30; + REGISTER Simd result_31; + REGISTER Simd result_32; // 20 left - REGISTER vComplex Chi_00; // two spinor; 6 regs - REGISTER vComplex Chi_01; - REGISTER vComplex Chi_02; + REGISTER Simd Chi_00; // two spinor; 6 regs + REGISTER Simd Chi_01; + REGISTER Simd Chi_02; - REGISTER vComplex Chi_10; - REGISTER vComplex Chi_11; - REGISTER vComplex Chi_12; // 14 left + REGISTER Simd Chi_10; + REGISTER Simd Chi_11; + REGISTER Simd Chi_12; // 14 left - REGISTER vComplex UChi_00; // two spinor; 6 regs - REGISTER vComplex UChi_01; - REGISTER vComplex UChi_02; + REGISTER Simd UChi_00; // two spinor; 6 regs + REGISTER Simd UChi_01; + REGISTER Simd UChi_02; - REGISTER vComplex UChi_10; - REGISTER vComplex UChi_11; - REGISTER vComplex UChi_12; // 8 left + REGISTER Simd UChi_10; + REGISTER Simd UChi_11; + REGISTER Simd UChi_12; // 8 left - REGISTER vComplex U_00; // two rows of U matrix - REGISTER vComplex U_10; - REGISTER vComplex U_20; - REGISTER vComplex U_01; - REGISTER vComplex U_11; - REGISTER vComplex U_21; // 2 reg left. + REGISTER Simd U_00; // two rows of U matrix + REGISTER Simd U_10; + REGISTER Simd U_20; + REGISTER Simd U_01; + REGISTER Simd U_11; + REGISTER Simd U_21; // 2 reg left. #define Chimu_00 Chi_00 #define Chimu_01 Chi_01 @@ -360,11 +362,6 @@ void DiracOptHandDhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U, MULT_2SPIN(Xp); } XP_RECON; - // std::cout< > &buf, - int ss,int sU,const LatticeFermion &in, LatticeFermion &out) +template +void WilsonKernels >::DiracOptHandDhopSiteDag(CartesianStencil &st,DoubledGaugeField &U, + std::vector > &buf, + int ss,int sU,const FermionField &in, FermionField &out) { - REGISTER vComplex result_00; // 12 regs on knc - REGISTER vComplex result_01; - REGISTER vComplex result_02; + REGISTER Simd result_00; // 12 regs on knc + REGISTER Simd result_01; + REGISTER Simd result_02; - REGISTER vComplex result_10; - REGISTER vComplex result_11; - REGISTER vComplex result_12; + REGISTER Simd result_10; + REGISTER Simd result_11; + REGISTER Simd result_12; - REGISTER vComplex result_20; - REGISTER vComplex result_21; - REGISTER vComplex result_22; + REGISTER Simd result_20; + REGISTER Simd result_21; + REGISTER Simd result_22; - REGISTER vComplex result_30; - REGISTER vComplex result_31; - REGISTER vComplex result_32; // 20 left + REGISTER Simd result_30; + REGISTER Simd result_31; + REGISTER Simd result_32; // 20 left - REGISTER vComplex Chi_00; // two spinor; 6 regs - REGISTER vComplex Chi_01; - REGISTER vComplex Chi_02; + REGISTER Simd Chi_00; // two spinor; 6 regs + REGISTER Simd Chi_01; + REGISTER Simd Chi_02; - REGISTER vComplex Chi_10; - REGISTER vComplex Chi_11; - REGISTER vComplex Chi_12; // 14 left + REGISTER Simd Chi_10; + REGISTER Simd Chi_11; + REGISTER Simd Chi_12; // 14 left - REGISTER vComplex UChi_00; // two spinor; 6 regs - REGISTER vComplex UChi_01; - REGISTER vComplex UChi_02; + REGISTER Simd UChi_00; // two spinor; 6 regs + REGISTER Simd UChi_01; + REGISTER Simd UChi_02; - REGISTER vComplex UChi_10; - REGISTER vComplex UChi_11; - REGISTER vComplex UChi_12; // 8 left + REGISTER Simd UChi_10; + REGISTER Simd UChi_11; + REGISTER Simd UChi_12; // 8 left - REGISTER vComplex U_00; // two rows of U matrix - REGISTER vComplex U_10; - REGISTER vComplex U_20; - REGISTER vComplex U_01; - REGISTER vComplex U_11; - REGISTER vComplex U_21; // 2 reg left. + REGISTER Simd U_00; // two rows of U matrix + REGISTER Simd U_10; + REGISTER Simd U_20; + REGISTER Simd U_01; + REGISTER Simd U_11; + REGISTER Simd U_21; // 2 reg left. #define Chimu_00 Chi_00 #define Chimu_01 Chi_01 @@ -752,7 +744,7 @@ void DiracOptHandDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, TP_RECON_ACCUM; { - vSpinColourVector & ref (out._odata[ss]); + SiteSpinor & ref (out._odata[ss]); vstream(ref()(0)(0),result_00*(-0.5)); vstream(ref()(0)(1),result_01*(-0.5)); vstream(ref()(0)(2),result_02*(-0.5)); @@ -767,4 +759,5 @@ void DiracOptHandDhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U, vstream(ref()(3)(2),result_32*(-0.5)); } } +#endif }} diff --git a/lib/qcd/action/pseudofermion/TwoFlavour.h b/lib/qcd/action/pseudofermion/TwoFlavour.h index 22f761fe..635c1ca6 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavour.h +++ b/lib/qcd/action/pseudofermion/TwoFlavour.h @@ -95,12 +95,13 @@ namespace Grid{ //////////////////////////////////////////////////////////////////////// // Two flavour pseudofermion action for any dop //////////////////////////////////////////////////////////////////////// - template - class TwoFlavourPseudoFermionAction : public Action { + template + class TwoFlavourPseudoFermionAction : public Action { private: +#include - FermionOperator & FermOp;// the basic operator + FermionOperator & FermOp;// the basic operator OperatorFunction &DerivativeSolver; @@ -112,7 +113,7 @@ namespace Grid{ ///////////////////////////////////////////////// // Pass in required objects. ///////////////////////////////////////////////// - TwoFlavourPseudoFermionAction(FermionOperator &Op, + TwoFlavourPseudoFermionAction(FermionOperator &Op, OperatorFunction & DS, OperatorFunction & AS ) : FermOp(Op), DerivativeSolver(DS), ActionSolver(AS), Phi(Op.FermionGrid()) { @@ -158,7 +159,7 @@ namespace Grid{ FermionField X(FermOp.FermionGrid()); FermionField Y(FermOp.FermionGrid()); - MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); X=zero; ActionSolver(MdagMOp,Phi,X); MdagMOp.Op(X,Y); @@ -183,7 +184,7 @@ namespace Grid{ FermionField Y(FermOp.FermionGrid()); GaugeField tmp(FermOp.GaugeGrid()); - MdagMLinearOperator ,FermionField> MdagMOp(FermOp); + MdagMLinearOperator ,FermionField> MdagMOp(FermOp); X=zero; DerivativeSolver(MdagMOp,Phi,X); diff --git a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h index 63dec41a..35a208c8 100644 --- a/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h +++ b/lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h @@ -4,13 +4,18 @@ namespace Grid{ namespace QCD{ - template - class SchurDifferentiableOperator : public SchurDiagMooeeOperator + template + class SchurDifferentiableOperator : public SchurDiagMooeeOperator,typename Impl::FermionField> { + public: +#include + public: + typedef FermionOperator Matrix; + SchurDifferentiableOperator (Matrix &Mat) : SchurDiagMooeeOperator(Mat) {}; - void MpcDeriv(LatticeGaugeField &Force,const FermionField &U,const FermionField &V) { + void MpcDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { GridBase *fgrid = this->_Mat.FermionGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); @@ -28,8 +33,8 @@ namespace Grid{ assert(U.checkerboard==Odd); assert(V.checkerboard==V.checkerboard); - LatticeGaugeField ForceO(ucbgrid); - LatticeGaugeField ForceE(ucbgrid); + GaugeField ForceO(ucbgrid); + GaugeField ForceE(ucbgrid); // X^dag Der_oe MeeInv Meo Y // Use Mooee as nontrivial but gauge field indept @@ -48,7 +53,7 @@ namespace Grid{ } - void MpcDagDeriv(LatticeGaugeField &Force,const FermionField &U,const FermionField &V) { + void MpcDagDeriv(GaugeField &Force,const FermionField &U,const FermionField &V) { GridBase *fgrid = this->_Mat.FermionGrid(); GridBase *fcbgrid = this->_Mat.FermionRedBlackGrid(); @@ -66,8 +71,8 @@ namespace Grid{ assert(V.checkerboard==Odd); assert(V.checkerboard==V.checkerboard); - LatticeGaugeField ForceO(ucbgrid); - LatticeGaugeField ForceE(ucbgrid); + GaugeField ForceO(ucbgrid); + GaugeField ForceE(ucbgrid); // X^dag Der_oe MeeInv Meo Y // Use Mooee as nontrivial but gauge field indept @@ -91,12 +96,15 @@ namespace Grid{ //////////////////////////////////////////////////////////////////////// // Two flavour pseudofermion action for any EO prec dop //////////////////////////////////////////////////////////////////////// - template - class TwoFlavourEvenOddPseudoFermionAction : public Action { + template + class TwoFlavourEvenOddPseudoFermionAction : public Action { + + public: +#include private: - FermionOperator & FermOp;// the basic operator + FermionOperator & FermOp;// the basic operator OperatorFunction &DerivativeSolver; @@ -109,7 +117,7 @@ namespace Grid{ ///////////////////////////////////////////////// // Pass in required objects. ///////////////////////////////////////////////// - TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, + TwoFlavourEvenOddPseudoFermionAction(FermionOperator &Op, OperatorFunction & DS, OperatorFunction & AS ) : @@ -140,7 +148,7 @@ namespace Grid{ pickCheckerboard(Even,etaEven,eta); pickCheckerboard(Odd,etaOdd,eta); - SchurDifferentiableOperator,FermionField> PCop(FermOp); + SchurDifferentiableOperator PCop(FermOp); FermOp.ImportGauge(U); @@ -163,7 +171,7 @@ namespace Grid{ FermionField X(FermOp.FermionRedBlackGrid()); FermionField Y(FermOp.FermionRedBlackGrid()); - SchurDifferentiableOperator,FermionField> PCop(FermOp); + SchurDifferentiableOperator PCop(FermOp); X=zero; ActionSolver(PCop,PhiOdd,X); @@ -195,7 +203,7 @@ namespace Grid{ FermionField Y(FermOp.FermionRedBlackGrid()); GaugeField tmp(FermOp.GaugeGrid()); - SchurDifferentiableOperator,FermionField> PCop(FermOp); + SchurDifferentiableOperator PCop(FermOp); X=zero; DerivativeSolver(PCop,PhiOdd,X); diff --git a/lib/qcd/spin/Dirac.h b/lib/qcd/spin/Dirac.h index 3208d0b5..65d8d6b5 100644 --- a/lib/qcd/spin/Dirac.h +++ b/lib/qcd/spin/Dirac.h @@ -4,7 +4,6 @@ namespace Grid{ namespace QCD { - const int SpinorIndex = 2; class Gamma { diff --git a/lib/qcd/spin/TwoSpinor.h b/lib/qcd/spin/TwoSpinor.h index 919f0f76..94ba9978 100644 --- a/lib/qcd/spin/TwoSpinor.h +++ b/lib/qcd/spin/TwoSpinor.h @@ -42,18 +42,14 @@ namespace QCD { * -i 0 0 0 */ - - template strong_inline void - spProjXp (iVector &hspin,const iVector &fspin) + // To fail is not to err (Cryptic clue: suggest to Google SFINAE ;) ) + template > = 0> strong_inline void spProjXp (iVector &hspin,const iVector &fspin) { - // To fail is not to err (Cryptic clue: suggest to Google SFINAE ;) ) - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)+timesI(fspin(3)); hspin(1)=fspin(1)+timesI(fspin(2)); } - template strong_inline void spProjXm (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjXm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)-timesI(fspin(3)); hspin(1)=fspin(1)-timesI(fspin(2)); } @@ -62,15 +58,13 @@ namespace QCD { // 0 0 1 0 [1] +- [2] // 0 1 0 0 // -1 0 0 0 - template strong_inline void spProjYp (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjYp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)-fspin(3); hspin(1)=fspin(1)+fspin(2); } - template strong_inline void spProjYm (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjYm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)+fspin(3); hspin(1)=fspin(1)-fspin(2); } @@ -80,15 +74,14 @@ namespace QCD { * -i 0 0 0 * 0 i 0 0 */ - template strong_inline void spProjZp (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjZp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)+timesI(fspin(2)); hspin(1)=fspin(1)-timesI(fspin(3)); } - template strong_inline void spProjZm (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjZm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)-timesI(fspin(2)); hspin(1)=fspin(1)+timesI(fspin(3)); } @@ -98,15 +91,15 @@ namespace QCD { * 1 0 0 0 * 0 1 0 0 */ - template strong_inline void spProjTp (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjTp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)+fspin(2); hspin(1)=fspin(1)+fspin(3); } - template strong_inline void spProjTm (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProjTm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0)-fspin(2); hspin(1)=fspin(1)-fspin(3); } @@ -117,32 +110,33 @@ namespace QCD { * 0 0 0 -1 */ - template strong_inline void spProj5p (iVector &hspin,const iVector &fspin) + template > = 0> strong_inline void spProj5p (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(0); hspin(1)=fspin(1); } - template strong_inline void spProj5m (iVector &hspin,const iVector &fspin) + + template > = 0> strong_inline void spProj5m (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; hspin(0)=fspin(2); hspin(1)=fspin(3); } // template strong_inline void fspProj5p (iVector &rfspin,const iVector &fspin) - template strong_inline void spProj5p (iVector &rfspin,const iVector &fspin) + template > = 0> strong_inline void spProj5p (iVector &rfspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; rfspin(0)=fspin(0); rfspin(1)=fspin(1); rfspin(2)=zero; rfspin(3)=zero; } // template strong_inline void fspProj5m (iVector &rfspin,const iVector &fspin) - template strong_inline void spProj5m (iVector &rfspin,const iVector &fspin) + template > = 0> strong_inline void spProj5m (iVector &rfspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; rfspin(0)=zero; rfspin(1)=zero; rfspin(2)=fspin(2); @@ -158,33 +152,33 @@ namespace QCD { * 0 -i 0 0 -i[1]+-[2] == -i ([0]+-i[3]) = -i (1) * -i 0 0 0 -i[0]+-[3] == -i ([1]+-i[2]) = -i (0) */ - template strong_inline void spReconXp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconXp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=timesMinusI(hspin(1)); fspin(3)=timesMinusI(hspin(0)); } - template strong_inline void spReconXm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconXm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=timesI(hspin(1)); fspin(3)=timesI(hspin(0)); } - template strong_inline void accumReconXp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconXp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)-=timesI(hspin(1)); fspin(3)-=timesI(hspin(0)); } - template strong_inline void accumReconXm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconXm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)+=timesI(hspin(1)); @@ -196,33 +190,33 @@ namespace QCD { // 0 1 0 0 == 1(1) // -1 0 0 0 ==-1(0) - template strong_inline void spReconYp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconYp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)= hspin(1); fspin(3)=-hspin(0);//Unary minus? } - template strong_inline void spReconYm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconYm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=-hspin(1); fspin(3)= hspin(0); } - template strong_inline void accumReconYp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconYp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)+=hspin(1); fspin(3)-=hspin(0); } - template strong_inline void accumReconYm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconYm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)-=hspin(1); @@ -235,33 +229,33 @@ namespace QCD { * -i 0 0 0 => -i (0) * 0 i 0 0 => i (1) */ - template strong_inline void spReconZp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconZp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=timesMinusI(hspin(0)); fspin(3)=timesI(hspin(1)); } - template strong_inline void spReconZm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconZm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)= timesI(hspin(0)); fspin(3)=timesMinusI(hspin(1)); } - template strong_inline void accumReconZp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconZp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)-=timesI(hspin(0)); fspin(3)+=timesI(hspin(1)); } - template strong_inline void accumReconZm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconZm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)+=timesI(hspin(0)); @@ -273,33 +267,33 @@ namespace QCD { * 1 0 0 0 => (0) * 0 1 0 0 => (1) */ - template strong_inline void spReconTp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconTp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=hspin(0); fspin(3)=hspin(1); } - template strong_inline void spReconTm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spReconTm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0); fspin(1)=hspin(1); fspin(2)=-hspin(0); fspin(3)=-hspin(1); } - template strong_inline void accumReconTp (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconTp (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)+=hspin(0); fspin(3)+=hspin(1); } - template strong_inline void accumReconTm (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumReconTm (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0); fspin(1)+=hspin(1); fspin(2)-=hspin(0); @@ -311,31 +305,31 @@ namespace QCD { * 0 0 -1 0 * 0 0 0 -1 */ - template strong_inline void spRecon5p (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spRecon5p (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=hspin(0)+hspin(0); // add is lower latency than mul fspin(1)=hspin(1)+hspin(1); // probably no measurable diffence though fspin(2)=zero; fspin(3)=zero; } - template strong_inline void spRecon5m (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void spRecon5m (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)=zero; fspin(1)=zero; fspin(2)=hspin(0)+hspin(0); fspin(3)=hspin(1)+hspin(1); } - template strong_inline void accumRecon5p (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumRecon5p (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(0)+=hspin(0)+hspin(0); fspin(1)+=hspin(1)+hspin(1); } - template strong_inline void accumRecon5m (iVector &fspin,const iVector &hspin) + template > = 0> strong_inline void accumRecon5m (iVector &fspin,const iVector &hspin) { - typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; + //typename std::enable_if,SpinorIndex>::value,iVector >::type *SFINAE; fspin(2)+=hspin(0)+hspin(0); fspin(3)+=hspin(1)+hspin(1); } @@ -347,21 +341,19 @@ namespace QCD { ////////// // Xp ////////// - template strong_inline void spProjXp (iScalar &hspin,const iScalar &fspin) + template = 0> strong_inline void spProjXp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; - spProjXp(hspin._internal,fspin._internal); - } - template strong_inline void spProjXp (iVector &hspin,iVector &fspin) - { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjXp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjXp (iScalar &hspin,const iScalar &fspin) + { + spProjXp(hspin._internal,fspin._internal); + } + template strong_inline void spProjXp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconXp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconXp(hspin._internal,fspin._internal); } - template strong_inline void spReconXp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconXp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconXp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconXp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconXp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconXp(hspin._internal,fspin._internal); } - template strong_inline void accumReconXp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconXp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconXp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconXp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjXm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjXm(hspin._internal,fspin._internal); } - template strong_inline void spProjXm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjXm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjXm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjXm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconXm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconXm(hspin._internal,fspin._internal); } - template strong_inline void spReconXm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconXm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconXm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconXm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconXm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconXm(hspin._internal,fspin._internal); } - template strong_inline void accumReconXm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconXm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconXm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconXm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjYp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjYp(hspin._internal,fspin._internal); } - template strong_inline void spProjYp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjYp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjYp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjYp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconYp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconYp(hspin._internal,fspin._internal); } - template strong_inline void spReconYp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconYp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconYp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconYp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconYp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconYp(hspin._internal,fspin._internal); } - template strong_inline void accumReconYp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconYp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconYp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconYp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjYm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjYm(hspin._internal,fspin._internal); } - template strong_inline void spProjYm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjYm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjYm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjYm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconYm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconYm(hspin._internal,fspin._internal); } - template strong_inline void spReconYm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconYm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,const iVector >::type *temp; for(int i=0;i strong_inline void spReconYm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconYm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconYm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconYm(hspin._internal,fspin._internal); } - template strong_inline void accumReconYm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconYm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconYm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconYm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjZp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjZp(hspin._internal,fspin._internal); } - template strong_inline void spProjZp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjZp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjZp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjZp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconZp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconZp(hspin._internal,fspin._internal); } - template strong_inline void spReconZp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconZp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconZp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconZp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconZp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconZp(hspin._internal,fspin._internal); } - template strong_inline void accumReconZp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconZp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconZp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconZp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjZm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjZm(hspin._internal,fspin._internal); } - template strong_inline void spProjZm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjZm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjZm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjZm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconZm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconZm(hspin._internal,fspin._internal); } - template strong_inline void spReconZm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconZm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconZm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconZm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconZm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconZm(hspin._internal,fspin._internal); } - template strong_inline void accumReconZm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconZm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconZm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconZm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjTp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjTp(hspin._internal,fspin._internal); } - template strong_inline void spProjTp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjTp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjTp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjTp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconTp (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconTp(hspin._internal,fspin._internal); } - template strong_inline void spReconTp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconTp (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconTp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconTp (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconTp (iScalar &hspin,const iScalar &fspin) + template strong_inline void accumReconTp (iScalar &hspin, iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconTp(hspin._internal,fspin._internal); } - template strong_inline void accumReconTp (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconTp (iVector &hspin, const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconTp (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconTp (iMatrix &hspin, const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProjTm (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProjTm(hspin._internal,fspin._internal); } - template strong_inline void spProjTm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProjTm (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProjTm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProjTm (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spReconTm (iScalar &hspin,const iScalar &fspin) + template strong_inline void spReconTm (iScalar &hspin, const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spReconTm(hspin._internal,fspin._internal); } - template strong_inline void spReconTm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spReconTm (iVector &hspin, const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spReconTm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spReconTm (iMatrix &hspin, const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumReconTm (iScalar &hspin,const iScalar &fspin) + template strong_inline void accumReconTm (iScalar &hspin, const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumReconTm(hspin._internal,fspin._internal); } - template strong_inline void accumReconTm (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumReconTm (iVector &hspin, const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumReconTm (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumReconTm (iMatrix &hspin, const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProj5p (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProj5p(hspin._internal,fspin._internal); } - template strong_inline void spProj5p (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProj5p (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProj5p (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProj5p (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spRecon5p (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spRecon5p(hspin._internal,fspin._internal); } - template strong_inline void spRecon5p (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spRecon5p (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spRecon5p (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spRecon5p (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumRecon5p (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumRecon5p(hspin._internal,fspin._internal); } - template strong_inline void accumRecon5p (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumRecon5p (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumRecon5p (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumRecon5p (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void fspProj5p (iScalar &hspin,const iScalar &fspin) template strong_inline void spProj5p (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProj5p(hspin._internal,fspin._internal); } // template strong_inline void fspProj5p (iVector &hspin,iVector &fspin) - template strong_inline void spProj5p (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProj5p (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void fspProj5p (iMatrix &hspin,iMatrix &fspin) - template strong_inline void spProj5p (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProj5p (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spProj5m (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProj5m(hspin._internal,fspin._internal); } - template strong_inline void spProj5m (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProj5m (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spProj5m (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProj5m (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void spRecon5m (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spRecon5m(hspin._internal,fspin._internal); } - template strong_inline void spRecon5m (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spRecon5m (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void spRecon5m (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spRecon5m (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void accumRecon5m (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; accumRecon5m(hspin._internal,fspin._internal); } - template strong_inline void accumRecon5m (iVector &hspin,iVector &fspin) + template = 0> strong_inline void accumRecon5m (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void accumRecon5m (iMatrix &hspin,iMatrix &fspin) + template strong_inline void accumRecon5m (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i strong_inline void fspProj5m (iScalar &hspin,const iScalar &fspin) template strong_inline void spProj5m (iScalar &hspin,const iScalar &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iScalar >::type *temp; spProj5m(hspin._internal,fspin._internal); } // template strong_inline void fspProj5m (iVector &hspin,iVector &fspin) - template strong_inline void spProj5m (iVector &hspin,iVector &fspin) + template = 0> strong_inline void spProj5m (iVector &hspin,const iVector &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iVector >::type *temp; for(int i=0;i strong_inline void fspProj5m (iMatrix &hspin,iMatrix &fspin) - template strong_inline void spProj5m (iMatrix &hspin,iMatrix &fspin) + template strong_inline void spProj5m (iMatrix &hspin,const iMatrix &fspin) { - typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; + //typename std::enable_if,SpinorIndex>::notvalue,iMatrix >::type *temp; for(int i=0;i(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + TestCGinversions(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; std::cout<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + TestCGinversions(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0); + TestCGinversions(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0); + TestCGinversions(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + ShamirZolotarevFermionR Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); + TestCGinversions(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonCayleyTanhFermionR Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestCGinversions(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonCayleyZolotarevFermionR Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); + TestCGinversions(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); } diff --git a/tests/Test_cayley_coarsen_support.cc b/tests/Test_cayley_coarsen_support.cc index 3ddf6205..77e77307 100644 --- a/tests/Test_cayley_coarsen_support.cc +++ b/tests/Test_cayley_coarsen_support.cc @@ -67,8 +67,8 @@ int main (int argc, char ** argv) RealD mass=0.5; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); HermIndefOp.Op(src,ref); HermIndefOp.OpDiag(src,result); @@ -89,7 +89,7 @@ int main (int argc, char ** argv) std::cout< HermDefOp(Ddwf); + MdagMLinearOperator HermDefOp(Ddwf); typedef Aggregation Subspace; Subspace Aggregates(Coarse5d,FGrid); Aggregates.CreateSubspaceRandom(RNG5); diff --git a/tests/Test_cayley_even_odd.cc b/tests/Test_cayley_even_odd.cc index d0c7b74a..c36567cd 100644 --- a/tests/Test_cayley_even_odd.cc +++ b/tests/Test_cayley_even_odd.cc @@ -49,35 +49,35 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5 =1.8; std::cout<(Ddwf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + DomainWallFermionR 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; std::cout<(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c); + TestWhat(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0); + TestWhat(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + ScaledShamirFermionR Dsham(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,2.0); + TestWhat(Dsham,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + ShamirZolotarevFermionR Dshamz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); + TestWhat(Dshamz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonCayleyTanhFermionR Dov(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestWhat(Dov,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonCayleyZolotarevFermionR Dovz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,2.0); + TestWhat(Dovz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); } diff --git a/tests/Test_cayley_ldop_cr.cc b/tests/Test_cayley_ldop_cr.cc index 48d7c9dc..2414d6f1 100644 --- a/tests/Test_cayley_ldop_cr.cc +++ b/tests/Test_cayley_ldop_cr.cc @@ -55,8 +55,8 @@ int main (int argc, char ** argv) std::cout< HermIndefOp(Ddwf); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); const int nbasis = 8; @@ -67,7 +67,7 @@ int main (int argc, char ** argv) std::cout< HermDefOp(Ddwf); + MdagMLinearOperator HermDefOp(Ddwf); Subspace Aggregates(Coarse5d,FGrid); Aggregates.CreateSubspace(RNG5,HermDefOp); diff --git a/tests/Test_cf_coarsen_support.cc b/tests/Test_cf_coarsen_support.cc index caa862be..d5cb7dfc 100644 --- a/tests/Test_cf_coarsen_support.cc +++ b/tests/Test_cf_coarsen_support.cc @@ -48,8 +48,8 @@ int main (int argc, char ** argv) RealD M5=1.8; { - OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); - HermitianLinearOperator HermIndefOp(Dcf); + OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + HermitianLinearOperator HermIndefOp(Dcf); HermIndefOp.Op(src,ref); HermIndefOp.OpDiag(src,result); @@ -65,8 +65,8 @@ int main (int argc, char ** argv) } { - OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); - HermitianLinearOperator HermIndefOp(Dpf); + OverlapWilsonPartialFractionTanhFermionR Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + HermitianLinearOperator HermIndefOp(Dpf); HermIndefOp.Op(src,ref); HermIndefOp.OpDiag(src,result); diff --git a/tests/Test_cf_cr_unprec.cc b/tests/Test_cf_cr_unprec.cc index 90acb503..7cb238ea 100644 --- a/tests/Test_cf_cr_unprec.cc +++ b/tests/Test_cf_cr_unprec.cc @@ -44,14 +44,14 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.8; - OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); ConjugateResidual MCR(1.0e-8,10000); - MdagMLinearOperator HermPosDefOp(Dcf); + MdagMLinearOperator HermPosDefOp(Dcf); MCR(HermPosDefOp,src,result); - HermitianLinearOperator HermIndefOp(Dcf); + HermitianLinearOperator HermIndefOp(Dcf); MCR(HermIndefOp,src,result); Grid_finalize(); diff --git a/tests/Test_contfrac_cg.cc b/tests/Test_contfrac_cg.cc index da6a81b9..abf46754 100644 --- a/tests/Test_contfrac_cg.cc +++ b/tests/Test_contfrac_cg.cc @@ -75,21 +75,21 @@ int main (int argc, char ** argv) std::cout<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestCGinversions(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonContFracZolotarevFermionR Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0); + TestCGinversions(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonPartialFractionTanhFermionR Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestCGinversions(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonPartialFractionZolotarevFermionR Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0); + TestCGinversions(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); diff --git a/tests/Test_contfrac_even_odd.cc b/tests/Test_contfrac_even_odd.cc index c8e5b093..7fcfc0f4 100644 --- a/tests/Test_contfrac_even_odd.cc +++ b/tests/Test_contfrac_even_odd.cc @@ -50,20 +50,20 @@ int main (int argc, char ** argv) RealD M5 =1.8; std::cout<(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonContFracTanhFermionR Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestWhat(Dcf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonContFracZolotarevFermionR Dcfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0); + TestWhat(Dcfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonPartialFractionTanhFermionR Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + TestWhat(Dpf,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); std::cout<(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); + OverlapWilsonPartialFractionZolotarevFermionR Dpfz(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,0.1,6.0); + TestWhat(Dpfz,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5); Grid_finalize(); } diff --git a/tests/Test_contfrac_force.cc b/tests/Test_contfrac_force.cc index a39d596d..58b42693 100644 --- a/tests/Test_contfrac_force.cc +++ b/tests/Test_contfrac_force.cc @@ -42,7 +42,7 @@ int main (int argc, char ** argv) //////////////////////////////////// RealD mass=0.01; RealD M5=1.8; - OverlapWilsonContFracTanhFermion Dcf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + OverlapWilsonContFracTanhFermionR Dcf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); Dcf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p diff --git a/tests/Test_dwf_cg_prec.cc b/tests/Test_dwf_cg_prec.cc index 633d280f..a39b6529 100644 --- a/tests/Test_dwf_cg_prec.cc +++ b/tests/Test_dwf_cg_prec.cc @@ -45,14 +45,14 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); LatticeFermion src_o(FrbGrid); LatticeFermion result_o(FrbGrid); pickCheckerboard(Odd,src_o,src); result_o=zero; - SchurDiagMooeeOperator HermOpEO(Ddwf); + SchurDiagMooeeOperator HermOpEO(Ddwf); ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,result_o); diff --git a/tests/Test_dwf_cg_schur.cc b/tests/Test_dwf_cg_schur.cc index ac3989d9..e9db375e 100644 --- a/tests/Test_dwf_cg_schur.cc +++ b/tests/Test_dwf_cg_schur.cc @@ -43,7 +43,7 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); ConjugateGradient CG(1.0e-8,10000); SchurRedBlackDiagMooeeSolve SchurSolver(CG); diff --git a/tests/Test_dwf_cg_unprec.cc b/tests/Test_dwf_cg_unprec.cc index 055e1bef..55d0c8d0 100644 --- a/tests/Test_dwf_cg_unprec.cc +++ b/tests/Test_dwf_cg_unprec.cc @@ -43,9 +43,9 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - MdagMLinearOperator HermOp(Ddwf); + MdagMLinearOperator HermOp(Ddwf); ConjugateGradient CG(1.0e-8,10000); CG(HermOp,src,result); diff --git a/tests/Test_dwf_cr_unprec.cc b/tests/Test_dwf_cr_unprec.cc index c2b75853..2d791ce2 100644 --- a/tests/Test_dwf_cr_unprec.cc +++ b/tests/Test_dwf_cr_unprec.cc @@ -50,12 +50,12 @@ int main (int argc, char ** argv) RealD mass=0.5; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - MdagMLinearOperator HermOp(Ddwf); + MdagMLinearOperator HermOp(Ddwf); MCR(HermOp,src,result); - Gamma5R5HermitianLinearOperator g5HermOp(Ddwf); + Gamma5R5HermitianLinearOperator g5HermOp(Ddwf); MCR(g5HermOp,src,result); diff --git a/tests/Test_dwf_even_odd.cc b/tests/Test_dwf_even_odd.cc index 4e3d9799..96853b4c 100644 --- a/tests/Test_dwf_even_odd.cc +++ b/tests/Test_dwf_even_odd.cc @@ -58,7 +58,7 @@ int main (int argc, char ** argv) RealD mass=0.1; RealD M5 =1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); LatticeFermion src_e (FrbGrid); LatticeFermion src_o (FrbGrid); @@ -187,7 +187,7 @@ int main (int argc, char ** argv) RealD t1,t2; - SchurDiagMooeeOperator HermOpEO(Ddwf); + SchurDiagMooeeOperator HermOpEO(Ddwf); HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); diff --git a/tests/Test_dwf_force.cc b/tests/Test_dwf_force.cc index a3dbb0f2..ac6fd26a 100644 --- a/tests/Test_dwf_force.cc +++ b/tests/Test_dwf_force.cc @@ -42,7 +42,7 @@ int main (int argc, char ** argv) //////////////////////////////////// RealD mass=0.01; RealD M5=1.8; - DomainWallFermion Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); Ddwf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p diff --git a/tests/Test_dwf_fpgcr.cc b/tests/Test_dwf_fpgcr.cc index 6c589c6c..6c9645b9 100644 --- a/tests/Test_dwf_fpgcr.cc +++ b/tests/Test_dwf_fpgcr.cc @@ -52,19 +52,19 @@ int main (int argc, char ** argv) RealD mass=0.5; RealD M5=1.8; - DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); std::cout< HermOp(Ddwf); + MdagMLinearOperator HermOp(Ddwf); result=zero; PGCR(HermOp,src,result); std::cout< g5HermOp(Ddwf); + Gamma5R5HermitianLinearOperator g5HermOp(Ddwf); result=zero; PGCR(g5HermOp,src,result); diff --git a/tests/Test_dwf_hdcr.cc b/tests/Test_dwf_hdcr.cc index 6ff26425..45647936 100644 --- a/tests/Test_dwf_hdcr.cc +++ b/tests/Test_dwf_hdcr.cc @@ -401,7 +401,7 @@ int main (int argc, char ** argv) std::cout< HermDefOp(Ddwf); + MdagMLinearOperator HermDefOp(Ddwf); Subspace Aggregates(Coarse5d,FGrid); // Aggregates.CreateSubspace(RNG5,HermDefOp,nbasis); assert ( (nbasis & 0x1)==0); @@ -437,7 +437,7 @@ int main (int argc, char ** argv) std::cout< HermIndefOp(Ddwf); + Gamma5R5HermitianLinearOperator HermIndefOp(Ddwf); CoarsenedMatrix LDOp(*Coarse5d); LDOp.CoarsenOperator(FGrid,HermIndefOp,Aggregates); @@ -467,7 +467,7 @@ int main (int argc, char ** argv) std::cout< Precon(Aggregates, LDOp,HermIndefOp,Ddwf); + MultiGridPreconditioner Precon(Aggregates, LDOp,HermIndefOp,Ddwf); TrivialPrecon simple; std::cout< HermOpEO(Ddwf); + SchurDiagMooeeOperator HermOpEO(Ddwf); ConjugateGradient pCG(1.0e-8,10000); LatticeFermion src_o(FrbGrid); diff --git a/tests/Test_hmc_EOWilsonFermionGauge.cc b/tests/Test_hmc_EOWilsonFermionGauge.cc index 1c90933e..66ec118b 100644 --- a/tests/Test_hmc_EOWilsonFermionGauge.cc +++ b/tests/Test_hmc_EOWilsonFermionGauge.cc @@ -25,13 +25,12 @@ int main (int argc, char ** argv) WilsonGaugeAction Waction(5.6); Real mass=-0.77; - WilsonFermion FermOp(U,Fine,RBFine,mass); + WilsonFermionR FermOp(U,Fine,RBFine,mass); ConjugateGradient CG(1.0e-8,10000); ConjugateGradient CGmd(1.0e-6,10000); - TwoFlavourEvenOddPseudoFermionAction - WilsonNf2(FermOp,CG,CG); + TwoFlavourEvenOddPseudoFermionAction WilsonNf2(FermOp,CG,CG); //Collect actions ActionLevel Level1; diff --git a/tests/Test_hmc_WilsonFermionGauge.cc b/tests/Test_hmc_WilsonFermionGauge.cc index e3bbf393..ef598062 100644 --- a/tests/Test_hmc_WilsonFermionGauge.cc +++ b/tests/Test_hmc_WilsonFermionGauge.cc @@ -25,13 +25,12 @@ int main (int argc, char ** argv) WilsonGaugeAction Waction(5.6); Real mass=-0.77; - WilsonFermion FermOp(U,Fine,RBFine,mass); + WilsonFermionR FermOp(U,Fine,RBFine,mass); ConjugateGradient CG(1.0e-8,10000); ConjugateGradient CGmd(1.0e-6,10000); - TwoFlavourPseudoFermionAction - WilsonNf2(FermOp,CG,CG); + TwoFlavourPseudoFermionAction WilsonNf2(FermOp,CG,CG); //Collect actions ActionLevel Level1; diff --git a/tests/Test_partfrac_force.cc b/tests/Test_partfrac_force.cc index 05925e55..6e543d3c 100644 --- a/tests/Test_partfrac_force.cc +++ b/tests/Test_partfrac_force.cc @@ -42,7 +42,7 @@ int main (int argc, char ** argv) //////////////////////////////////// RealD mass=0.01; RealD M5=1.8; - OverlapWilsonPartialFractionTanhFermion Dpf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); + OverlapWilsonPartialFractionTanhFermionR Dpf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0); Dpf.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p diff --git a/tests/Test_wilson_cg_prec.cc b/tests/Test_wilson_cg_prec.cc index c59699e4..756f8b6f 100644 --- a/tests/Test_wilson_cg_prec.cc +++ b/tests/Test_wilson_cg_prec.cc @@ -42,7 +42,7 @@ int main (int argc, char ** argv) } RealD mass=0.5; - WilsonFermion Dw(Umu,Grid,RBGrid,mass); + WilsonFermionR Dw(Umu,Grid,RBGrid,mass); // HermitianOperator HermOp(Dw); // ConjugateGradient CG(1.0e-8,10000); @@ -53,7 +53,7 @@ int main (int argc, char ** argv) pickCheckerboard(Odd,src_o,src); result_o=zero; - SchurDiagMooeeOperator HermOpEO(Dw); + SchurDiagMooeeOperator HermOpEO(Dw); ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o,result_o); diff --git a/tests/Test_wilson_cg_schur.cc b/tests/Test_wilson_cg_schur.cc index 1484a6a8..2174d644 100644 --- a/tests/Test_wilson_cg_schur.cc +++ b/tests/Test_wilson_cg_schur.cc @@ -37,7 +37,7 @@ int main (int argc, char ** argv) LatticeFermion resid(&Grid); RealD mass=0.5; - WilsonFermion Dw(Umu,Grid,RBGrid,mass); + WilsonFermionR Dw(Umu,Grid,RBGrid,mass); ConjugateGradient CG(1.0e-8,10000); SchurRedBlackDiagMooeeSolve SchurSolver(CG); diff --git a/tests/Test_wilson_cg_unprec.cc b/tests/Test_wilson_cg_unprec.cc index 5984cac9..9269aeaa 100644 --- a/tests/Test_wilson_cg_unprec.cc +++ b/tests/Test_wilson_cg_unprec.cc @@ -40,9 +40,9 @@ int main (int argc, char ** argv) } RealD mass=0.5; - WilsonFermion Dw(Umu,Grid,RBGrid,mass); + WilsonFermionR Dw(Umu,Grid,RBGrid,mass); - MdagMLinearOperator HermOp(Dw); + MdagMLinearOperator HermOp(Dw); ConjugateGradient CG(1.0e-8,10000); CG(HermOp,src,result); diff --git a/tests/Test_wilson_cr_unprec.cc b/tests/Test_wilson_cr_unprec.cc index 468f0d5a..c7e0bc33 100644 --- a/tests/Test_wilson_cr_unprec.cc +++ b/tests/Test_wilson_cr_unprec.cc @@ -43,9 +43,9 @@ int main (int argc, char ** argv) } RealD mass=0.5; - WilsonFermion Dw(Umu,Grid,RBGrid,mass); + WilsonFermionR Dw(Umu,Grid,RBGrid,mass); - MdagMLinearOperator HermOp(Dw); + MdagMLinearOperator HermOp(Dw); ConjugateResidual MCR(1.0e-8,10000); diff --git a/tests/Test_wilson_even_odd.cc b/tests/Test_wilson_even_odd.cc index 66ac89dc..37ba40f1 100644 --- a/tests/Test_wilson_even_odd.cc +++ b/tests/Test_wilson_even_odd.cc @@ -62,7 +62,7 @@ int main (int argc, char ** argv) RealD mass=0.1; - WilsonFermion Dw(Umu,Grid,RBGrid,mass); + WilsonFermionR Dw(Umu,Grid,RBGrid,mass); LatticeFermion src_e (&RBGrid); LatticeFermion src_o (&RBGrid); @@ -179,7 +179,7 @@ int main (int argc, char ** argv) pickCheckerboard(Odd ,phi_o,phi); RealD t1,t2; - SchurDiagMooeeOperator HermOpEO(Dw); + SchurDiagMooeeOperator HermOpEO(Dw); HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2); HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2); diff --git a/tests/Test_wilson_force.cc b/tests/Test_wilson_force.cc index b08ccd45..1b22c2e4 100644 --- a/tests/Test_wilson_force.cc +++ b/tests/Test_wilson_force.cc @@ -38,7 +38,7 @@ int main (int argc, char ** argv) // Unmodified matrix element //////////////////////////////////// RealD mass=-4.0; //kills the diagonal term - WilsonFermion Dw (U, Grid,RBGrid,mass); + WilsonFermionR Dw (U, Grid,RBGrid,mass); Dw.M (phi,Mphi); ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p @@ -88,7 +88,7 @@ int main (int argc, char ** argv) } std::cout << GridLogMessage <<"Initial mom hamiltonian is "<< Hmom <