portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2024-09-20 01:05:38 +01:00
Code Releases Activity

Domain wall fermions now invert ; have the basis set up for

Browse Source 
Tanh/Zolo * (Cayley/PartFrac/ContFrac) * (Mobius/Shamir/Wilson)
Approx        Representation               Kernel.

All are done with space-time taking part in checkerboarding, Ls uncheckerboarded

Have only so far tested the Domain Wall limit of mobius, and at that only checked
that it
i)  Inverts
ii) 5dim DW == Ls copies of 4dim D2
iii) MeeInv Mee == 1
iv) Meo+Mee+Moe+Moo == M unprec.
v) MpcDagMpc is hermitan
vi) Mdag is the adjoint of M between stochastic vectors.

That said, the RB schur solve, RB MpcDagMpc solve, Unprec solve
all converge and the true residual becomes small; so pretty good tests.
This commit is contained in:
Peter Boyle 2015-06-02 16:57:12 +01:00
parent c851d0e705
commit 3845f267cb
36 changed files with 1500 additions and 220 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
benchmarks/Grid_dwf.cc
View File

@ -24,43 +24,28 @@ int main (int argc, char ** argv)
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> latt4 = GridDefaultLatt();
std::vector<int> simd4 = GridDefaultSimd(Nd,vComplexF::Nsimd());
std::vector<int> mpi4 = GridDefaultMpi();
assert(latt4.size()==4 );
assert(simd4.size()==4 );
assert(mpi4.size() ==4 );
const int Ls=1;
std::vector<int> latt5({Ls,latt4[0],latt4[1],latt4[2],latt4[3]});
std::vector<int> simd5({1 ,simd4[0],simd4[1],simd4[2],simd4[3]});
std::vector<int> mpi5({1 , mpi4[0], mpi4[1], mpi4[2], mpi4[3]});
std::vector<int> cb5({0,1,1,1,1}); // Checkerboard 4d only
int cbd=1; // use dim-1 to reduce
// Four dim grid for gauge field U
GridCartesian UGrid(latt4,simd4,mpi4);
GridRedBlackCartesian UrbGrid(&UGrid);
// Five dim grid for fermions F
GridCartesian FGrid(latt5,simd5,mpi5);
GridRedBlackCartesian FrbGrid(latt5,simd5,mpi5,cb5,cbd);
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(&FGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeFermion src (&FGrid); random(RNG5,src);
LatticeFermion result(&FGrid); result=zero;
LatticeFermion ref(&FGrid); ref=zero;
LatticeFermion tmp(&FGrid);
LatticeFermion err(&FGrid);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeFermion src (FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
ColourMatrix cm = Complex(1.0,0.0);
GridParallelRNG RNG4(&UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(&UGrid); random(RNG4,Umu);
LatticeGaugeField Umu5d(&FGrid);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu5d(FGrid);
// replicate across fifth dimension
for(int ss=0;ss<Umu._grid->oSites();ss++){
@ -72,7 +57,7 @@ int main (int argc, char ** argv)
////////////////////////////////////
// Naive wilson implementation
////////////////////////////////////
std::vector<LatticeColourMatrix> U(4,&FGrid);
std::vector<LatticeColourMatrix> U(4,FGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = peekIndex<LorentzIndex>(Umu5d,mu);
}
@ -93,17 +78,17 @@ int main (int argc, char ** argv)
}
RealD mass=0.1;
FiveDimWilsonFermion Dw(Umu,FGrid,FrbGrid,UGrid,UrbGrid,mass);
RealD M5 =1.8;
DomainWallFermion Dw(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
std::cout << "Calling Dw"<<std::endl;
int ncall=1000;
int ncall=10;
double t0=usecond();
for(int i=0;i<ncall;i++){
Dw.Dhop(src,result,0);
}
double t1=usecond();
double volume=Ls; for(int mu=0;mu<Nd;mu++) volume=volume*latt4[mu];
double flops=1344*volume*ncall;
@ -141,5 +126,31 @@ int main (int argc, char ** argv)
err = ref-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
LatticeFermion src_e (FrbGrid);
LatticeFermion src_o (FrbGrid);
LatticeFermion r_e (FrbGrid);
LatticeFermion r_o (FrbGrid);
LatticeFermion r_eo (FGrid);
std::cout << "Calling Deo and Doe"<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Dw.DhopEO(src_o,r_e,DaggerNo);
Dw.DhopOE(src_e,r_o,DaggerNo);
Dw.Dhop(src,result,DaggerNo);
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
err = r_eo-result;
std::cout << "norm diff "<< norm2(err)<<std::endl;
pickCheckerboard(Even,src_e,err);
pickCheckerboard(Odd,src_o,err);
std::cout << "norm diff even "<< norm2(src_e)<<std::endl;
std::cout << "norm diff odd "<< norm2(src_o)<<std::endl;
Grid_finalize();
}

58
benchmarks/Grid_dwf_cg_prec.cc Normal file
View File

@ -0,0 +1,58 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template<class d>
struct scal {
d internal;
};
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
}
RealD mass=0.1;
RealD M5=1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd,src_o,src);
result_o=zero;
HermitianCheckerBoardedOperator<DomainWallFermion,LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOpEO,src_o,result_o);
Grid_finalize();
}

53
benchmarks/Grid_dwf_cg_schur.cc Normal file
View File

@ -0,0 +1,53 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template<class d>
struct scal {
d internal;
};
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
}
RealD mass=0.1;
RealD M5=1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
SchurRedBlackSolve<LatticeFermion> SchurSolver(CG);
SchurSolver(Ddwf,src,result);
Grid_finalize();
}

53
benchmarks/Grid_dwf_cg_unprec.cc Normal file
View File

@ -0,0 +1,53 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template<class d>
struct scal {
d internal;
};
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
}
RealD mass=0.1;
RealD M5=1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
HermitianOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
CG(HermOp,src,result);
Grid_finalize();
}

207
benchmarks/Grid_dwf_even_odd.cc Normal file
View File

@ -0,0 +1,207 @@
#include <Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template<class d>
struct scal {
d internal;
};
Gamma::GammaMatrix Gmu [] = {
Gamma::GammaX,
Gamma::GammaY,
Gamma::GammaZ,
Gamma::GammaT
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeFermion src (FGrid); random(RNG5,src);
LatticeFermion phi (FGrid); random(RNG5,phi);
LatticeFermion chi (FGrid); random(RNG5,chi);
LatticeFermion result(FGrid); result=zero;
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid); tmp=zero;
LatticeFermion err(FGrid); tmp=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)
Umu=zero;
for(int nn=0;nn<Nd;nn++){
random(RNG4,U[nn]);
if ( nn>0 )
U[nn]=zero;
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
RealD mass=0.1;
RealD M5 =1.8;
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
LatticeFermion src_e (FrbGrid);
LatticeFermion src_o (FrbGrid);
LatticeFermion r_e (FrbGrid);
LatticeFermion r_o (FrbGrid);
LatticeFermion r_eo (FGrid);
LatticeFermion r_eeoo(FGrid);
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ddwf.Meooe(src_e,r_o); std::cout<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<"Applied Moe"<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
Ddwf.Mooee(src_e,r_e); std::cout<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<"Applied Moo"<<std::endl;
setCheckerboard(r_eeoo,r_e);
setCheckerboard(r_eeoo,r_o);
r_eo=r_eo+r_eeoo;
Ddwf.M(src,ref);
// std::cout << r_eo<<std::endl;
// std::cout << ref <<std::endl;
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(FGrid);
cerr = localInnerProduct(err,err);
// std::cout << cerr<<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
LatticeFermion chi_e (FrbGrid);
LatticeFermion chi_o (FrbGrid);
LatticeFermion dchi_e (FrbGrid);
LatticeFermion dchi_o (FrbGrid);
LatticeFermion phi_e (FrbGrid);
LatticeFermion phi_o (FrbGrid);
LatticeFermion dphi_e (FrbGrid);
LatticeFermion dphi_o (FrbGrid);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
Ddwf.Meooe(chi_e,dchi_o);
Ddwf.Meooe(chi_o,dchi_e);
Ddwf.MeooeDag(phi_e,dphi_o);
Ddwf.MeooeDag(phi_o,dphi_e);
ComplexD pDce = innerProduct(phi_e,dchi_e);
ComplexD pDco = innerProduct(phi_o,dchi_o);
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ddwf.Mooee(chi_e,src_e);
Ddwf.MooeeInv(src_e,phi_e);
Ddwf.Mooee(chi_o,src_o);
Ddwf.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ddwf.MooeeDag(chi_e,src_e);
Ddwf.MooeeInvDag(src_e,phi_e);
Ddwf.MooeeDag(chi_o,src_o);
Ddwf.MooeeInvDag(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout << "norm diff "<< norm2(err)<< std::endl;
std::cout<<"=============================================================="<<std::endl;
std::cout<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<"=============================================================="<<std::endl;
random(RNG5,phi);
random(RNG5,chi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
Ddwf.MpcDagMpc(chi_e,dchi_e,t1,t2);
Ddwf.MpcDagMpc(chi_o,dchi_o,t1,t2);
Ddwf.MpcDagMpc(phi_e,dphi_e,t1,t2);
Ddwf.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}

5
benchmarks/Grid_wilson_evenodd.cc
View File

@ -68,8 +68,6 @@ int main (int argc, char ** argv)
LatticeFermion r_o (&RBGrid);
LatticeFermion r_eo (&Grid);
const int Even=0;
const int Odd=1;
std::cout<<"=========================================================="<<std::endl;
std::cout<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
std::cout<<"=========================================================="<<std::endl;
@ -79,12 +77,11 @@ int main (int argc, char ** argv)
Dw.Meooe(src_e,r_o); std::cout<<"Applied Meo"<<std::endl;
Dw.Meooe(src_o,r_e); std::cout<<"Applied Moe"<<std::endl;
Dw.Dhop (src,ref,0);
Dw.Dhop (src,ref,DaggerNo);
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
ref = (-0.5)*ref;
err= ref - r_eo;
std::cout << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;

46
benchmarks/Makefile.am
View File

@ -5,14 +5,34 @@ AM_LDFLAGS = -L$(top_builddir)/lib
#
# Test code
#
bin_PROGRAMS = Grid_wilson Grid_comms Grid_memory_bandwidth Grid_su3 Grid_wilson_cg_unprec Grid_wilson_evenodd Grid_wilson_cg_prec Grid_wilson_cg_schur Grid_dwf
bin_PROGRAMS = \
Grid_comms \
Grid_memory_bandwidth \
Grid_su3 \
Grid_wilson \
Grid_wilson_evenodd \
Grid_wilson_cg_unprec \
Grid_wilson_cg_prec \
Grid_wilson_cg_schur \
Grid_dwf\
Grid_dwf_even_odd\
Grid_dwf_cg_unprec\
Grid_dwf_cg_prec\
Grid_dwf_cg_schur
Grid_comms_SOURCES = Grid_comms.cc
Grid_comms_LDADD = -lGrid
Grid_su3_SOURCES = Grid_su3.cc Grid_su3_test.cc Grid_su3_expr.cc
Grid_su3_LDADD = -lGrid
Grid_memory_bandwidth_SOURCES = Grid_memory_bandwidth.cc
Grid_memory_bandwidth_LDADD = -lGrid
Grid_wilson_SOURCES = Grid_wilson.cc
Grid_wilson_LDADD = -lGrid
Grid_dwf_SOURCES = Grid_dwf.cc
Grid_dwf_LDADD = -lGrid
Grid_wilson_evenodd_SOURCES = Grid_wilson_evenodd.cc
Grid_wilson_evenodd_LDADD = -lGrid
@ -25,12 +45,18 @@ Grid_wilson_cg_prec_LDADD = -lGrid
Grid_wilson_cg_schur_SOURCES = Grid_wilson_cg_schur.cc
Grid_wilson_cg_schur_LDADD = -lGrid
Grid_comms_SOURCES = Grid_comms.cc
Grid_comms_LDADD = -lGrid
Grid_dwf_SOURCES = Grid_dwf.cc
Grid_dwf_LDADD = -lGrid
Grid_su3_SOURCES = Grid_su3.cc Grid_su3_test.cc Grid_su3_expr.cc
Grid_su3_LDADD = -lGrid
Grid_dwf_even_odd_SOURCES = Grid_dwf_even_odd.cc
Grid_dwf_even_odd_LDADD = -lGrid
Grid_memory_bandwidth_SOURCES = Grid_memory_bandwidth.cc
Grid_memory_bandwidth_LDADD = -lGrid
Grid_dwf_cg_unprec_SOURCES = Grid_dwf_cg_unprec.cc
Grid_dwf_cg_unprec_LDADD = -lGrid
Grid_dwf_cg_prec_SOURCES = Grid_dwf_cg_prec.cc
Grid_dwf_cg_prec_LDADD = -lGrid
Grid_dwf_cg_schur_SOURCES = Grid_dwf_cg_schur.cc
Grid_dwf_cg_schur_LDADD = -lGrid

2
lib/Grid_init.cc
View File

@ -143,7 +143,7 @@ void Grid_init(int *argc,char ***argv)
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--dslash-opt") ){
WilsonFermion::HandOptDslash=1;
FiveDimWilsonFermion::HandOptDslash=1;
WilsonFermion5D::HandOptDslash=1;
}
if( GridCmdOptionExists(*argv,*argv+*argc,"--lebesgue") ){
LebesgueOrder::UseLebesgueOrder=1;

16
lib/Makefile.am
View File

@ -20,13 +20,17 @@ libGrid_a_SOURCES = \
stencil/Grid_stencil_common.cc \
algorithms/approx/Zolotarev.cc \
algorithms/approx/Remez.cc \
qcd/action/fermion/FiveDimWilsonFermion.cc\
qcd/action/fermion/WilsonFermion.cc\
qcd/SpaceTimeGrid.cc\
qcd/Dirac.cc\
qcd/action/fermion/WilsonKernels.cc\
qcd/action/fermion/WilsonKernelsHand.cc\
qcd/Dirac.cc\
qcd/action/fermion/WilsonFermion.cc\
qcd/action/fermion/WilsonFermion5D.cc\
qcd/action/fermion/CayleyFermion5D.cc \
qcd/action/fermion/ContinuedFractionFermion5D.cc \
$(extra_sources)
# qcd/action/fermion/PartialFractionFermion5D.cc \
#
# Include files
#
@ -95,11 +99,11 @@ nobase_include_HEADERS=\
./math/Grid_math_transpose.h\
./parallelIO/GridNerscIO.h\
./qcd/action/Actions.h\
./qcd/action/fermion/FermionAction.h\
./qcd/action/fermion/FiveDimWilsonFermion.h\
./qcd/action/fermion/FermionOperator.h\
./qcd/action/fermion/WilsonCompressor.h\
./qcd/action/fermion/WilsonFermion.h\
./qcd/action/fermion/WilsonKernels.h\
./qcd/action/fermion/WilsonFermion.h\
./qcd/action/fermion/WilsonFermion5D.h\
./qcd/Dirac.h\
./qcd/QCD.h\
./qcd/TwoSpinor.h\

32
lib/algorithms/LinearOperator.h
View File

@ -125,39 +125,7 @@ namespace Grid {
};
*/
// Chroma interface defining GaugeAction
/*
template<typename P, typename Q> class GaugeAction
virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
virtual GaugeState<P,Q>* createState(const Q& q) const
virtual const GaugeBC<P,Q>& getGaugeBC() const
virtual const Set& getSet(void) const = 0;
virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const
virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
typedef multi1d<LatticeColorMatrix> P;
typedef multi1d<LatticeColorMatrix> Q;
virtual void staple(LatticeColorMatrix& result,
const Handle< GaugeState<P,Q> >& state,
int mu, int cb) const = 0;
*/
// Chroma interface defining FermionAction
/*
template<typename T, typename P, typename Q> class FermAct4D : public FermionAction<T,P,Q>
virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
class DiffFermAct4D : public FermAct4D<T,P,Q>
virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
*/
}
#endif

9
lib/algorithms/approx/Zolotarev.cc
View File

@ -58,6 +58,8 @@
/* Compute the partial fraction expansion coefficients (alpha) from the
* factored form */
namespace Grid {
namespace Approx {
static void construct_partfrac(izd *z) {
int dn = z -> dn, dd = z -> dd, type = z -> type;
@ -291,7 +293,7 @@ static void sncndnFK(INTERNAL_PRECISION u, INTERNAL_PRECISION k,
* Set type = 0 for the Zolotarev approximation, which is zero at x = 0, and
* type = 1 for the approximation which is infinite at x = 0. */
zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
zolotarev_data* grid_zolotarev(PRECISION epsilon, int n, int type) {
INTERNAL_PRECISION A, c, cp, kp, ksq, sn, cn, dn, Kp, Kj, z, z0, t, M, F,
l, invlambda, xi, xisq, *tv, s, opl;
int m, czero, ts;
@ -412,7 +414,7 @@ zolotarev_data* bfm_zolotarev(PRECISION epsilon, int n, int type) {
return zd;
}
zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
zolotarev_data* grid_higham(PRECISION epsilon, int n) {
INTERNAL_PRECISION A, M, c, cp, z, z0, t, epssq;
int m, czero;
zolotarev_data *zd;
@ -502,6 +504,7 @@ zolotarev_data* bfm_higham(PRECISION epsilon, int n) {
free(d);
return zd;
}
}}
#ifdef TEST
@ -707,4 +710,6 @@ int main(int argc, char** argv) {
return EXIT_SUCCESS;
}
#endif /* TEST */

9
lib/algorithms/approx/Zolotarev.h
View File

@ -1,7 +1,8 @@
/* -*- Mode: C; comment-column: 22; fill-column: 79; -*- */
#ifdef __cplusplus
extern "C" {
namespace Grid {
namespace Approx {
#endif
#define HVERSION Header Time-stamp: <14-OCT-2004 09:26:51.00 adk@MISSCONTRARY>
@ -76,10 +77,10 @@ typedef struct {
* zolotarev_data structure. The arguments must satisfy the constraints that
* epsilon > 0, n > 0, and type = 0 or 1. */
ZOLOTAREV_DATA* bfm_higham(PRECISION epsilon, int n) ;
ZOLOTAREV_DATA* bfm_zolotarev(PRECISION epsilon, int n, int type);
ZOLOTAREV_DATA* grid_higham(PRECISION epsilon, int n) ;
ZOLOTAREV_DATA* grid_zolotarev(PRECISION epsilon, int n, int type);
#endif
#ifdef __cplusplus
}
}}
#endif

2
lib/cartesian/Grid_cartesian_base.h
View File

@ -21,7 +21,7 @@ public:
// Give Lattice access
template<class object> friend class Lattice;
GridBase(std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
GridBase(const std::vector<int> & processor_grid) : CartesianCommunicator(processor_grid) {};
// Physics Grid information.

6
lib/cartesian/Grid_cartesian_full.h
View File

@ -27,9 +27,9 @@ public:
virtual int CheckerBoardShift(int source_cb,int dim,int shift, int osite){
return shift;
}
GridCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid
GridCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid
) : GridBase(processor_grid)
{
///////////////////////

24
lib/cartesian/Grid_cartesian_red_black.h
View File

@ -81,28 +81,28 @@ public:
}
};
GridRedBlackCartesian(GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {};
GridRedBlackCartesian(const GridBase *base) : GridRedBlackCartesian(base->_fdimensions,base->_simd_layout,base->_processors) {};
GridRedBlackCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid,
std::vector<int> &checker_dim_mask,
GridRedBlackCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid,
const std::vector<int> &checker_dim_mask,
int checker_dim
) : GridBase(processor_grid)
{
Init(dimensions,simd_layout,processor_grid,checker_dim_mask,checker_dim);
}
GridRedBlackCartesian(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid) : GridBase(processor_grid)
GridRedBlackCartesian(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid) : GridBase(processor_grid)
{
std::vector<int> checker_dim_mask(dimensions.size(),1);
Init(dimensions,simd_layout,processor_grid,checker_dim_mask,0);
}
void Init(std::vector<int> &dimensions,
std::vector<int> &simd_layout,
std::vector<int> &processor_grid,
std::vector<int> &checker_dim_mask,
void Init(const std::vector<int> &dimensions,
const std::vector<int> &simd_layout,
const std::vector<int> &processor_grid,
const std::vector<int> &checker_dim_mask,
int checker_dim)
{
///////////////////////

2
lib/communicator/Grid_communicator_base.h
View File

@ -27,7 +27,7 @@ class CartesianCommunicator {
#endif
// Constructor
CartesianCommunicator(std::vector<int> &pdimensions_in);
CartesianCommunicator(const std::vector<int> &pdimensions_in);
// Wraps MPI_Cart routines
void ShiftedRanks(int dim,int shift,int & source, int & dest);

2
lib/communicator/Grid_communicator_mpi.cc
View File

@ -5,7 +5,7 @@ namespace Grid {
// Should error check all MPI calls.
CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
{
_ndimension = processors.size();
std::vector<int> periodic(_ndimension,1);

2
lib/communicator/Grid_communicator_none.cc
View File

@ -1,7 +1,7 @@
#include "Grid.h"
namespace Grid {
CartesianCommunicator::CartesianCommunicator(std::vector<int> &processors)
CartesianCommunicator::CartesianCommunicator(const std::vector<int> &processors)
{
_processors = processors;
_ndimension = processors.size();

113
lib/qcd/LinalgUtils.h Normal file
View File

@ -0,0 +1,113 @@
#ifndef GRID_QCD_LINALG_UTILS_H
#define GRID_QCD_LINALG_UTILS_H
namespace Grid{
namespace QCD{
////////////////////////////////////////////////////////////////////////
//This file brings additional linear combination assist that is helpful
//to QCD such as chiral projectors and spin matrices applied to one of the inputs.
//These routines support five-D chiral fermions and contain s-subslice indexing
//on the 5d (rb4d) checkerboarded lattices
////////////////////////////////////////////////////////////////////////
template<class vobj>
void axpby_ssp(Lattice<vobj> &z, RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp = a*x._odata[ss+s]+b*y._odata[ss+sp];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void ag5xpby_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
multGamma5(tmp(),a*x._odata[ss+s]());
tmp = tmp + b*y._odata[ss+sp];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void axpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
multGamma5(tmp(),b*y._odata[ss+sp]());
tmp = tmp + a*x._odata[ss+s];
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void ag5xpbg5y_ssp(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp1;
vobj tmp2;
tmp1 = a*x._odata[ss+s]+b*y._odata[ss+sp];
multGamma5(tmp2(),tmp1());
vstream(z._odata[ss+s],tmp2);
}
}
template<class vobj>
void axpby_ssp_pminus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
spProj5m(tmp,y._odata[ss+sp]);
tmp = a*x._odata[ss+s]+b*tmp;
vstream(z._odata[ss+s],tmp);
}
}
template<class vobj>
void axpby_ssp_pplus(Lattice<vobj> &z,RealD a,const Lattice<vobj> &x,RealD b,const Lattice<vobj> &y,int s,int sp)
{
z.checkerboard = x.checkerboard;
conformable(x,y);
conformable(x,z);
GridBase *grid=x._grid;
int Ls = grid->_rdimensions[0];
PARALLEL_FOR_LOOP
for(int ss=0;ss<grid->oSites();ss+=Ls){ // adds Ls
vobj tmp;
spProj5p(tmp,y._odata[ss+sp]);
tmp = a*x._odata[ss+s]+b*tmp;
vstream(z._odata[ss+s],tmp);
}
}
}}
#endif

2
lib/qcd/QCD.h
View File

@ -307,8 +307,10 @@ namespace QCD {
} //namespace QCD
} // Grid
#include <qcd/SpaceTimeGrid.h>
#include <qcd/Dirac.h>
#include <qcd/TwoSpinor.h>
#include <qcd/LinalgUtils.h>
#include <qcd/action/Actions.h>
#endif

52
lib/qcd/SpaceTimeGrid.cc Normal file
View File

@ -0,0 +1,52 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
/////////////////////////////////////////////////////////////////
// Public interface
/////////////////////////////////////////////////////////////////
GridCartesian *SpaceTimeGrid::makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi)
{
return new GridCartesian(latt,simd,mpi);
}
GridRedBlackCartesian *SpaceTimeGrid::makeFourDimRedBlackGrid(const GridCartesian *FourDimGrid)
{
return new GridRedBlackCartesian(FourDimGrid);
}
GridCartesian *SpaceTimeGrid::makeFiveDimGrid(int Ls,const GridCartesian *FourDimGrid)
{
int N4=FourDimGrid->_ndimension;
std::vector<int> latt5(1,Ls);
std::vector<int> simd5(1,1);
std::vector<int> mpi5(1,1);
for(int d=0;d<N4;d++){
latt5.push_back(FourDimGrid->_fdimensions[d]);
simd5.push_back(FourDimGrid->_simd_layout[d]);
mpi5.push_back(FourDimGrid->_processors[d]);
}
return new GridCartesian(latt5,simd5,mpi5);
}
GridRedBlackCartesian *SpaceTimeGrid::makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid)
{
int N4=FourDimGrid->_ndimension;
int cbd=1;
std::vector<int> latt5(1,Ls);
std::vector<int> simd5(1,1);
std::vector<int> mpi5(1,1);
std::vector<int> cb5(1,0);
for(int d=0;d<N4;d++){
latt5.push_back(FourDimGrid->_fdimensions[d]);
simd5.push_back(FourDimGrid->_simd_layout[d]);
mpi5.push_back(FourDimGrid->_processors[d]);
cb5.push_back( 1);
}
return new GridRedBlackCartesian(latt5,simd5,mpi5,cb5,cbd);
}
}}

18
lib/qcd/SpaceTimeGrid.h Normal file
View File

@ -0,0 +1,18 @@
#ifndef GRID_QCD_SPACE_TIME_GRID_H
#define GRID_QCD_SPACE_TIME_GRID_H
namespace Grid {
namespace QCD {
class SpaceTimeGrid {
public:
static GridCartesian *makeFourDimGrid(const std::vector<int> & latt,const std::vector<int> &simd,const std::vector<int> &mpi);
static GridRedBlackCartesian *makeFourDimRedBlackGrid (const GridCartesian *FourDimGrid);
static GridCartesian *makeFiveDimGrid (int Ls,const GridCartesian *FourDimGrid);
static GridRedBlackCartesian *makeFiveDimRedBlackGrid(int Ls,const GridCartesian *FourDimGrid);
};
}}
#endif

78
lib/qcd/action/Actions.h
View File

@ -1,10 +1,80 @@
#ifndef GRID_QCD_ACTIONS_H
#define GRID_QCD_ACTIONS_H
#include <qcd/action/fermion/FermionAction.h>
#include <qcd/action/fermion/WilsonCompressor.h>
#include <qcd/action/fermion/WilsonKernels.h>
// Some reorganisation likely required as both Chroma and IroIro
// are separating the concept of the operator from that of action.
//
// The FermAction contains methods to create
//
// * Linear operators (Hermitian and non-hermitian) .. my LinearOperator
// * System solvers (Hermitian and non-hermitian) .. my OperatorFunction
// * MultiShift System solvers (Hermitian and non-hermitian) .. my OperatorFunction
////////////////////////////////////////////
// Abstract base interface
////////////////////////////////////////////
#include <qcd/action/fermion/FermionOperator.h>
////////////////////////////////////////////
// Utility functions
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonCompressor.h> //used by all wilson type fermions
#include <qcd/action/fermion/WilsonKernels.h> //used by all wilson type fermions
////////////////////////////////////////////
// 4D formulations
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonFermion.h>
#include <qcd/action/fermion/FiveDimWilsonFermion.h>
//#include <qcd/action/fermion/CloverFermion.h>
////////////////////////////////////////////
// 5D formulations
////////////////////////////////////////////
#include <qcd/action/fermion/WilsonFermion5D.h> // used by all 5d overlap types
#include <qcd/action/fermion/CayleyFermion5D.h>
#include <qcd/action/fermion/ContinuedFractionFermion5D.h>
//#include <qcd/action/fermion/PartialFraction.h>
#include <qcd/action/fermion/DomainWallFermion.h>
//#include <qcd/action/fermion/ScaledShamirCayleyTanh.h>
// Chroma interface defining FermionAction
/*
template<typename T, typename P, typename Q> class FermAct4D : public FermionAction<T,P,Q>
virtual LinearOperator<T>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinearOperator<T>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
virtual LinOpSystemSolver<T>* invLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMSystemSolver<T>* invMdagM(Handle< FermState<T,P,Q> > state,
virtual LinOpMultiSystemSolver<T>* mInvLinOp(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolver<T>* mInvMdagM(Handle< FermState<T,P,Q> > state,
virtual MdagMMultiSystemSolverAccumulate<T>* mInvMdagMAcc(Handle< FermState<T,P,Q> > state,
virtual SystemSolver<T>* qprop(Handle< FermState<T,P,Q> > state,
class DiffFermAct4D : public FermAct4D<T,P,Q>
virtual DiffLinearOperator<T,Q,P>* linOp(Handle< FermState<T,P,Q> > state) const = 0;
virtual DiffLinearOperator<T,Q,P>* lMdagM(Handle< FermState<T,P,Q> > state) const = 0;
*/
// Chroma interface defining GaugeAction
/*
template<typename P, typename Q> class GaugeAction
virtual const CreateGaugeState<P,Q>& getCreateState() const = 0;
virtual GaugeState<P,Q>* createState(const Q& q) const
virtual const GaugeBC<P,Q>& getGaugeBC() const
virtual const Set& getSet(void) const = 0;
virtual void deriv(P& result, const Handle< GaugeState<P,Q> >& state) const
virtual Double S(const Handle< GaugeState<P,Q> >& state) const = 0;
class LinearGaugeAction : public GaugeAction< multi1d<LatticeColorMatrix>, multi1d<LatticeColorMatrix> >
typedef multi1d<LatticeColorMatrix> P;
typedef multi1d<LatticeColorMatrix> Q;
virtual void staple(LatticeColorMatrix& result,
const Handle< GaugeState<P,Q> >& state,
int mu, int cb) const = 0;
*/
#endif

235
lib/qcd/action/fermion/CayleyFermion5D.cc Normal file
View File

@ -0,0 +1,235 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
CayleyFermion5D::CayleyFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5) :
WilsonFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_M5),
mass(_mass)
{
std::cout << "Constructing a CayleyFermion5D"<<std::endl;
}
// override multiply
RealD CayleyFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion Din(psi._grid);
// Assemble Din
for(int s=0;s<Ls;s++){
if ( s==0 ) {
// Din = bs psi[s] + cs[s] psi[s+1}
axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
// Din+= -mass*cs[s] psi[s+1}
axpby_ssp_pplus (Din,1.0,Din,-mass*cs[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(Din,bs[s],psi,-mass*cs[s],psi,s,0);
axpby_ssp_pplus (Din,1.0,Din,cs[s],psi,s,s-1);
} else {
axpby_ssp_pminus(Din,bs[s],psi,cs[s],psi,s,s+1);
axpby_ssp_pplus(Din,1.0,Din,cs[s],psi,s,s-1);
}
}
DW(Din,chi,DaggerNo);
// ((b D_W + D_w hop terms +1) on s-diag
axpby(chi,1.0,1.0,chi,psi);
for(int s=0;s<Ls;s++){
if ( s==0 ){
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,0);
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
} else {
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s-1);
}
}
return norm2(chi);
}
RealD CayleyFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &chi)
{
// Under adjoint
//D1+ D1- P- -> D1+^dag P+ D2-^dag
//D2- P+ D2+ P-D1-^dag D2+dag
LatticeFermion Din(psi._grid);
// Apply Dw
DW(psi,Din,DaggerYes);
for(int s=0;s<Ls;s++){
// Collect the terms in DW
// Chi = bs Din[s] + cs[s] Din[s+1}
// Chi+= -mass*cs[s] psi[s+1}
if ( s==0 ) {
axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,-mass*cs[Ls-1],Din,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus (chi,bs[s],Din,-mass*cs[0],Din,s,0);
axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
} else {
axpby_ssp_pplus (chi,bs[s],Din,cs[s+1],Din,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,cs[s-1],Din,s,s-1);
}
// Collect the terms indept of DW
if ( s==0 ){
axpby_ssp_pplus (chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,mass,psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus (chi,1.0,chi,mass,psi,s,0);
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
} else {
axpby_ssp_pplus(chi,1.0,chi,-1.0,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,-1.0,psi,s,s-1);
}
}
// ((b D_W + D_w hop terms +1) on s-diag
axpby (chi,1.0,1.0,chi,psi);
return norm2(chi);
}
// half checkerboard operations
void CayleyFermion5D::Meooe (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion tmp(psi._grid);
// Assemble the 5d matrix
for(int s=0;s<Ls;s++){
if ( s==0 ) {
// tmp = bs psi[s] + cs[s] psi[s+1}
// tmp+= -mass*cs[s] psi[s+1}
axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi ,s, s+1);
axpby_ssp_pplus(tmp,1.0,tmp,mass*ceo[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(tmp,beo[s],psi,mass*ceo[s],psi,s,0);
axpby_ssp_pplus(tmp,1.0,tmp,-ceo[s],psi,s,s-1);
} else {
axpby_ssp_pminus(tmp,beo[s],psi,-ceo[s],psi,s,s+1);
axpby_ssp_pplus (tmp,1.0,tmp,-ceo[s],psi,s,s-1);
}
}
// Apply 4d dslash
if ( psi.checkerboard == Odd ) {
DhopEO(tmp,chi,DaggerNo);
} else {
DhopOE(tmp,chi,DaggerNo);
}
}
void CayleyFermion5D::MeooeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion tmp(psi._grid);
// Apply 4d dslash
if ( psi.checkerboard == Odd ) {
DhopEO(psi,tmp,DaggerYes);
} else {
DhopOE(psi,tmp,DaggerYes);
}
// Assemble the 5d matrix
for(int s=0;s<Ls;s++){
if ( s==0 ) {
axpby_ssp_pplus(chi,beo[s],tmp, -ceo[s+1] ,tmp,s,s+1);
axpby_ssp_pminus(chi, 1.0,chi,mass*ceo[Ls-1],tmp,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus(chi,beo[s],tmp,mass*ceo[0],tmp,s,0);
axpby_ssp_pminus(chi,1.0,chi,-ceo[s-1],tmp,s,s-1);
} else {
axpby_ssp_pplus(chi,beo[s],tmp,-ceo[s+1],tmp,s,s+1);
axpby_ssp_pminus(chi,1.0 ,chi,-ceo[s-1],tmp,s,s-1);
}
}
}
void CayleyFermion5D::Mooee (const LatticeFermion &psi, LatticeFermion &chi)
{
for (int s=0;s<Ls;s++){
if ( s==0 ) {
axpby_ssp_pminus(chi,bee[s],psi ,-cee[s],psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,mass*cee[s],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pminus(chi,bee[s],psi,mass*cee[s],psi,s,0);
axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
} else {
axpby_ssp_pminus(chi,bee[s],psi,-cee[s],psi,s,s+1);
axpby_ssp_pplus (chi,1.0,chi,-cee[s],psi,s,s-1);
}
}
}
void CayleyFermion5D::MooeeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
for (int s=0;s<Ls;s++){
// Assemble the 5d matrix
if ( s==0 ) {
axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1] ,psi,s,s+1);
axpby_ssp_pminus(chi,1.0,chi,mass*cee[Ls-1],psi,s,Ls-1);
} else if ( s==(Ls-1)) {
axpby_ssp_pplus(chi,bee[s],psi,mass*cee[0],psi,s,0);
axpby_ssp_pminus(chi,1.0,chi,-cee[s-1],psi,s,s-1);
} else {
axpby_ssp_pplus(chi,bee[s],psi,-cee[s+1],psi,s,s+1);
axpby_ssp_pminus(chi,1.0 ,chi,-cee[s-1],psi,s,s-1);
}
}
}
void CayleyFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply (L^{\prime})^{-1}
axpby_ssp (chi,1.0,psi, 0.0,psi,0,0); // chi[0]=psi[0]
for (int s=1;s<Ls;s++){
axpby_ssp_pplus(chi,1.0,psi,-lee[s-1],chi,s,s-1);// recursion Psi[s] -lee P_+ chi[s-1]
}
// L_m^{-1}
for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
axpby_ssp_pminus(chi,1.0,chi,-leem[s],chi,Ls-1,s);
}
// U_m^{-1} D^{-1}
for (int s=0;s<Ls-1;s++){
// Chi[s] + 1/d chi[s]
axpby_ssp_pplus(chi,1.0/dee[s],chi,-ueem[s]/dee[Ls-1],chi,s,Ls-1);
}
axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable
// Apply U^{-1}
for (int s=Ls-2;s>=0;s--){
axpby_ssp_pminus (chi,1.0,chi,-uee[s],chi,s,s+1); // chi[Ls]
}
}
void CayleyFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply (U^{\prime})^{-dagger}
axpby_ssp (chi,1.0,psi, 0.0,psi,0,0); // chi[0]=psi[0]
for (int s=1;s<Ls;s++){
axpby_ssp_pminus(chi,1.0,psi,-uee[s-1],chi,s,s-1);
}
// U_m^{-\dagger}
for (int s=0;s<Ls-1;s++){
axpby_ssp_pplus(chi,1.0,chi,-ueem[s],chi,Ls-1,s);
}
// L_m^{-\dagger} D^{-dagger}
for (int s=0;s<Ls-1;s++){
axpby_ssp_pminus(chi,1.0/dee[s],chi,-leem[s]/dee[Ls-1],chi,s,Ls-1);
}
axpby_ssp(chi,1.0/dee[Ls-1],chi,0.0,chi,Ls-1,Ls-1); // Modest avoidable
// Apply L^{-dagger}
for (int s=Ls-2;s>=0;s--){
axpby_ssp_pplus (chi,1.0,chi,-lee[s],chi,s,s+1); // chi[Ls]
}
}
}
}

61
lib/qcd/action/fermion/CayleyFermion5D.h Normal file
View File

@ -0,0 +1,61 @@
#ifndef GRID_QCD_CAYLEY_FERMION_H
#define GRID_QCD_CAYLEY_FERMION_H
namespace Grid {
namespace QCD {
class CayleyFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// half checkerboard operations
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);
// protected:
Approx::zolotarev_data *zdata;
RealD mass;
// Cayley form Moebius (tanh and zolotarev)
std::vector<RealD> omega;
std::vector<RealD> bs; // S dependent coeffs
std::vector<RealD> cs;
std::vector<RealD> as;
// For preconditioning Cayley form
std::vector<RealD> bee;
std::vector<RealD> cee;
std::vector<RealD> aee;
std::vector<RealD> beo;
std::vector<RealD> ceo;
std::vector<RealD> aeo;
// LDU factorisation of the eeoo matrix
std::vector<RealD> lee;
std::vector<RealD> leem;
std::vector<RealD> uee;
std::vector<RealD> ueem;
std::vector<RealD> dee;
// Constructors
CayleyFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5);
};
}
}
#endif

119
lib/qcd/action/fermion/ContinuedFractionFermion5D.cc Normal file
View File

@ -0,0 +1,119 @@
#include <Grid.h>
namespace Grid {
namespace QCD {
RealD ContinuedFractionFermion5D::M (const LatticeFermion &psi, LatticeFermion &chi)
{
LatticeFermion D(psi._grid);
DW(psi,D,DaggerNo);
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==0 ) {
ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*scale,D,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
} else if ( s==(Ls-1) ){
RealD R=(1.0+mass)/(1.0-mass);
ag5xpby_ssp(chi,Beta[s]*scale,D,sqrt_cc[s-1],psi,s,s-1);
ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*scale,D,sqrt_cc[s],psi,s,s+1);
axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
}
sign=-sign;
}
return norm2(chi);
}
RealD ContinuedFractionFermion5D::Mdag (const LatticeFermion &psi, LatticeFermion &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::Meooe (const LatticeFermion &psi, LatticeFermion &chi)
{
Dhop(psi,chi,DaggerNo); // Dslash on diagonal. g5 Dslash is hermitian
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==(Ls-1) ){
ag5xpby_ssp(chi,Beta[s]*scale,chi,0.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*scale,chi,0.0,chi,s,s);
}
sign=-sign;
}
}
void ContinuedFractionFermion5D::MeooeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
Meooe(psi,chi);
}
void ContinuedFractionFermion5D::Mooee (const LatticeFermion &psi, LatticeFermion &chi)
{
double dw_diag = (4.0-this->M5)*scale;
int sign=1;
for(int s=0;s<Ls;s++){
if ( s==0 ) {
ag5xpby_ssp(chi,cc[0]*Beta[0]*sign*dw_diag,psi,sqrt_cc[0],psi,s,s+1); // Multiplies Dw by G5 so Hw
} else if ( s==(Ls-1) ){
// Drop the CC here.
double R=(1+this->mass)/(1-this->mass);
ag5xpby_ssp(chi,Beta[s]*dw_diag,psi,sqrt_cc[s-1],psi,s,s-1);
ag5xpby_ssp(chi,R,psi,1.0,chi,s,s);
} else {
ag5xpby_ssp(chi,cc[s]*Beta[s]*sign*dw_diag,psi,sqrt_cc[s],psi,s,s+1);
axpby_ssp(chi,1.0,chi,sqrt_cc[s-1],psi,s,s-1);
}
sign=-sign;
}
}
void ContinuedFractionFermion5D::MooeeDag (const LatticeFermion &psi, LatticeFermion &chi)
{
Mooee(psi,chi);
}
void ContinuedFractionFermion5D::MooeeInv (const LatticeFermion &psi, LatticeFermion &chi)
{
// Apply Linv
axpby_ssp(chi,1.0/cc_d[0],psi,0.0,psi,0,0);
for(int s=1;s<Ls;s++){
axpbg5y_ssp(chi,1.0/cc_d[s],psi,-1.0/See[s-1],chi,s,s-1);
}
// Apply Dinv
for(int s=0;s<Ls;s++){
ag5xpby_ssp(chi,1.0/See[s],chi,0.0,chi,s,s); //only appearance of See[0]
}
// Apply Uinv = (Linv)^T
axpby_ssp(chi,1.0/cc_d[Ls-1],chi,0.0,chi,this->Ls-1,this->Ls-1);
for(int s=Ls-2;s>=0;s--){
axpbg5y_ssp(chi,1.0/cc_d[s],chi,-1.0*cc_d[s+1]/See[s]/cc_d[s],chi,s,s+1);
}
}
void ContinuedFractionFermion5D::MooeeInvDag (const LatticeFermion &psi, LatticeFermion &chi)
{
MooeeInv(psi,chi);
}
// Constructors
ContinuedFractionFermion5D::ContinuedFractionFermion5D(
LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5) :
WilsonFermion5D(_Umu,
FiveDimGrid, FiveDimRedBlackGrid,
FourDimGrid, FourDimRedBlackGrid,M5),
mass(_mass)
{
}
}
}

53
lib/qcd/action/fermion/ContinuedFractionFermion5D.h Normal file
View File

@ -0,0 +1,53 @@
#ifndef GRID_QCD_CONTINUED_FRACTION_H
#define GRID_QCD_CONTINUED_FRACTION_H
namespace Grid {
namespace QCD {
class ContinuedFractionFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &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);
private:
Approx::zolotarev_data *zdata;
// Cont frac
RealD mass;
RealD R;
RealD scale;
std::vector<double> Beta;
std::vector<double> cc;;
std::vector<double> cc_d;;
std::vector<double> sqrt_cc;
std::vector<double> See;
std::vector<double> Aee;
// Constructors
ContinuedFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5);
};
}
}
#endif

118
lib/qcd/action/fermion/DomainWallFermion.h Normal file
View File

@ -0,0 +1,118 @@
#ifndef GRID_QCD_DOMAIN_WALL_FERMION_H
#define GRID_QCD_DOMAIN_WALL_FERMION_H
#include <Grid.h>
namespace Grid {
namespace QCD {
class DomainWallFermion : public CayleyFermion5D
{
public:
// Constructors
DomainWallFermion(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD _M5) :
CayleyFermion5D(_Umu,
FiveDimGrid,
FiveDimRedBlackGrid,
FourDimGrid,
FourDimRedBlackGrid,_mass,_M5)
{
RealD eps = 1.0;
zdata = Approx::grid_higham(eps,this->Ls);// eps is ignored for higham
assert(zdata->n==this->Ls);
///////////////////////////////////////////////////////////
// The Cayley coeffs (unprec)
///////////////////////////////////////////////////////////
this->omega.resize(this->Ls);
this->bs.resize(this->Ls);
this->cs.resize(this->Ls);
this->as.resize(this->Ls);
for(int i=0; i < this->Ls; i++){
this->as[i] = 1.0;
this->omega[i] = ((double)zdata -> gamma[i]);
double bb=1.0;
this->bs[i] = 0.5*(bb/(this->omega[i]) + 1.0);
this->cs[i] = 0.5*(bb/(this->omega[i]) - 1.0);
}
////////////////////////////////////////////////////////
// Constants for the preconditioned matrix Cayley form
////////////////////////////////////////////////////////
this->bee.resize(this->Ls);
this->cee.resize(this->Ls);
this->beo.resize(this->Ls);
this->ceo.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
this->bee[i]=as[i]*(bs[i]*(4.0-M5) +1.0);
this->cee[i]=as[i]*(1.0-cs[i]*(4.0-M5));
this->beo[i]=as[i]*bs[i];
this->ceo[i]=-as[i]*cs[i];
}
aee.resize(this->Ls);
aeo.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
aee[i]=cee[i];
aeo[i]=ceo[i];
}
//////////////////////////////////////////
// LDU decomposition of eeoo
//////////////////////////////////////////
dee.resize(this->Ls);
lee.resize(this->Ls);
leem.resize(this->Ls);
uee.resize(this->Ls);
ueem.resize(this->Ls);
for(int i=0;i<this->Ls;i++){
dee[i] = bee[i];
if ( i < this->Ls-1 ) {
lee[i] =-cee[i+1]/bee[i]; // sub-diag entry on the ith column
leem[i]=this->mass*cee[this->Ls-1]/bee[0];
for(int j=0;j<i;j++) leem[i]*= aee[j]/bee[j+1];
uee[i] =-aee[i]/bee[i]; // up-diag entry on the ith row
ueem[i]=this->mass;
for(int j=1;j<=i;j++) ueem[i]*= cee[j]/bee[j];
ueem[i]*= aee[0]/bee[0];
} else {
lee[i] =0.0;
leem[i]=0.0;
uee[i] =0.0;
ueem[i]=0.0;
}
}
{
double delta_d=mass*cee[this->Ls-1];
for(int j=0;j<this->Ls-1;j++) delta_d *= cee[j]/bee[j];
dee[this->Ls-1] += delta_d;
}
}
};
}
}
#endif

7
lib/qcd/action/fermion/FermionAction.h → lib/qcd/action/fermion/FermionOperator.h
View File

@ -1,5 +1,5 @@
#ifndef GRID_QCD_WILSON_DOP_H
#define GRID_QCD_WILSON_DOP_H
#ifndef GRID_QCD_FERMION_OPERATOR_H
#define GRID_QCD_FERMION_OPERATOR_H
namespace Grid {
@ -11,7 +11,7 @@ namespace Grid {
// Think about multiple representations
//////////////////////////////////////////////////////////////////////////////
template<class FermionField,class GaugeField>
class FermionAction : public CheckerBoardedSparseMatrixBase<FermionField>
class FermionOperator : public CheckerBoardedSparseMatrixBase<FermionField>
{
public:
@ -40,6 +40,7 @@ namespace Grid {
virtual void DhopOE(const FermionField &in, FermionField &out,int dag)=0;
virtual void DhopEO(const FermionField &in, FermionField &out,int dag)=0;
};
}

47
lib/qcd/action/fermion/PartialFractionFermion5D.cc Normal file
View File

@ -0,0 +1,47 @@
#ifndef GRID_QCD_PARTIAL_FRACTION_H
#define GRID_QCD_PARTIAL_FRACTION_H
namespace Grid {
namespace QCD {
class PartialFractionFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &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);
private:
zolotarev_data *zdata;
// Part frac
double R=(1+this->mass)/(1-this->mass);
std::vector<double> p;
std::vector<double> q;
// Constructors
PartialFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5);
};
}
}
#endif

49
lib/qcd/action/fermion/PartialFractionFermion5D.h Normal file
View File

@ -0,0 +1,49 @@
#ifndef GRID_QCD_PARTIAL_FRACTION_H
#define GRID_QCD_PARTIAL_FRACTION_H
namespace Grid {
namespace QCD {
class PartialFractionFermion5D : public WilsonFermion5D
{
public:
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &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);
private:
virtual void PartialFractionCoefficients(void);
zolotarev_data *zdata;
// Part frac
double R=(1+this->mass)/(1-this->mass);
std::vector<double> p;
std::vector<double> q;
// Constructors
PartialFractionFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
RealD _mass,RealD M5);
};
}
}
#endif

6
lib/qcd/action/fermion/WilsonFermion.cc
View File

@ -9,9 +9,9 @@ const std::vector<int> WilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
int WilsonFermion::HandOptDslash;
WilsonFermion::WilsonFermion(LatticeGaugeField &_Umu,
GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid,
double _mass) :
GridCartesian &Fgrid,
GridRedBlackCartesian &Hgrid,
RealD _mass) :
_grid(&Fgrid),
_cbgrid(&Hgrid),
Stencil (&Fgrid,npoint,Even,directions,displacements),

6
lib/qcd/action/fermion/WilsonFermion.h
View File

@ -5,7 +5,7 @@ namespace Grid {
namespace QCD {
class WilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
class WilsonFermion : public FermionOperator<LatticeFermion,LatticeGaugeField>
{
public:
@ -44,7 +44,7 @@ namespace Grid {
int dag);
// Constructor
WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass);
WilsonFermion(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,RealD _mass);
// DoubleStore
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@ -57,7 +57,7 @@ namespace Grid {
protected:
double mass;
RealD mass;
GridBase * _grid;
GridBase * _cbgrid;

106
lib/qcd/action/fermion/FiveDimWilsonFermion.cc → lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -4,18 +4,18 @@ namespace Grid {
namespace QCD {
// S-direction is INNERMOST and takes no part in the parity.
const std::vector<int> FiveDimWilsonFermion::directions ({1,2,3,4, 1, 2, 3, 4});
const std::vector<int> FiveDimWilsonFermion::displacements({1,1,1,1,-1,-1,-1,-1});
const std::vector<int> WilsonFermion5D::directions ({1,2,3,4, 1, 2, 3, 4});
const std::vector<int> WilsonFermion5D::displacements({1,1,1,1,-1,-1,-1,-1});
int FiveDimWilsonFermion::HandOptDslash;
int WilsonFermion5D::HandOptDslash;
// 5d lattice for DWF.
FiveDimWilsonFermion::FiveDimWilsonFermion(LatticeGaugeField &_Umu,
WilsonFermion5D::WilsonFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
double _mass) :
RealD _M5) :
_FiveDimGrid(&FiveDimGrid),
_FiveDimRedBlackGrid(&FiveDimRedBlackGrid),
_FourDimGrid(&FourDimGrid),
@ -23,7 +23,7 @@ namespace QCD {
Stencil (_FiveDimGrid,npoint,Even,directions,displacements),
StencilEven(_FiveDimRedBlackGrid,npoint,Even,directions,displacements), // source is Even
StencilOdd (_FiveDimRedBlackGrid,npoint,Odd ,directions,displacements), // source is Odd
mass(_mass),
M5(_M5),
Umu(_FourDimGrid),
UmuEven(_FourDimRedBlackGrid),
UmuOdd (_FourDimRedBlackGrid),
@ -70,7 +70,7 @@ namespace QCD {
pickCheckerboard(Even,UmuEven,Umu);
pickCheckerboard(Odd ,UmuOdd,Umu);
}
void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
void WilsonFermion5D::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
{
conformable(Uds._grid,GaugeGrid());
conformable(Umu._grid,GaugeGrid());
@ -82,60 +82,9 @@ void FiveDimWilsonFermion::DoubleStore(LatticeDoubledGaugeField &Uds,const Latti
pokeIndex<LorentzIndex>(Uds,U,mu+4);
}
}
RealD FiveDimWilsonFermion::M(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerNo);
return axpy_norm(out,5.0-M5,in,out);
}
RealD FiveDimWilsonFermion::Mdag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerYes);
return axpy_norm(out,5.0-M5,in,out);
}
void FiveDimWilsonFermion::Meooe(const LatticeFermion &in, LatticeFermion &out)
{
if ( in.checkerboard == Odd ) {
DhopEO(in,out,DaggerNo);
} else {
DhopOE(in,out,DaggerNo);
}
}
void FiveDimWilsonFermion::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
{
if ( in.checkerboard == Odd ) {
DhopEO(in,out,DaggerYes);
} else {
DhopOE(in,out,DaggerYes);
}
}
void FiveDimWilsonFermion::Mooee(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
out = (5.0-M5)*in;
return ;
}
void FiveDimWilsonFermion::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
Mooee(in,out);
}
void FiveDimWilsonFermion::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
out = (1.0/(5.0-M5))*in;
return ;
}
void FiveDimWilsonFermion::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
{
out.checkerboard = in.checkerboard;
MooeeInv(in,out);
}
void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopInternal(CartesianStencil & st, LebesgueOrder &lo,
LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
{
assert((dag==DaggerNo) ||(dag==DaggerYes));
@ -150,19 +99,21 @@ void FiveDimWilsonFermion::DhopInternal(CartesianStencil & st, LebesgueOrder &lo
// - 8 linear access unit stride streams per thread for Fermion for hw prefetchable.
if ( dag == DaggerYes ) {
if( HandOptDslash ) {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
//int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOptHand::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
}
}
} else {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOpt::DhopSiteDag(st,U,comm_buf,sF,sU,in,out);
}
@ -170,21 +121,22 @@ PARALLEL_FOR_LOOP
}
} else {
if( HandOptDslash ) {
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOptHand::DhopSite(st,U,comm_buf,sF,sU,in,out);
}
}
} else {
for(int ss=0;ss<U._grid->oSites();ss++){
int sU=lo.Reorder(ss);
PARALLEL_FOR_LOOP
for(int ss=0;ss<U._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
// int sU=lo.Reorder(ss);
int sU=ss;
int sF = s+Ls*sU;
DiracOpt::DhopSite(st,U,comm_buf,sF,sU,in,out);
}
@ -192,7 +144,7 @@ PARALLEL_FOR_LOOP
}
}
}
void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
conformable(in._grid,out._grid); // drops the cb check
@ -202,7 +154,7 @@ void FiveDimWilsonFermion::DhopOE(const LatticeFermion &in, LatticeFermion &out,
DhopInternal(StencilEven,LebesgueEvenOdd,UmuOdd,in,out,dag);
}
void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionRedBlackGrid()); // verifies half grid
conformable(in._grid,out._grid); // drops the cb check
@ -212,7 +164,7 @@ void FiveDimWilsonFermion::DhopEO(const LatticeFermion &in, LatticeFermion &out,
DhopInternal(StencilOdd,LebesgueEvenOdd,UmuEven,in,out,dag);
}
void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
void WilsonFermion5D::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag)
{
conformable(in._grid,FermionGrid()); // verifies full grid
conformable(in._grid,out._grid);
@ -221,8 +173,14 @@ void FiveDimWilsonFermion::Dhop(const LatticeFermion &in, LatticeFermion &out,in
DhopInternal(Stencil,Lebesgue,Umu,in,out,dag);
}
}}
void WilsonFermion5D::DW(const LatticeFermion &in, LatticeFermion &out,int dag)
{
out.checkerboard=in.checkerboard;
Dhop(in,out,dag); // -0.5 is included
axpy(out,4.0-M5,in,out);
}
}
}

31
lib/qcd/action/fermion/FiveDimWilsonFermion.h → lib/qcd/action/fermion/WilsonFermion5D.h
View File

@ -15,7 +15,7 @@ namespace Grid {
//
// [DIFFERS from original CPS red black implementation parity = (x+y+z+t+s)|2 ]
////////////////////////////////////////////////////////////////////////////////
class FiveDimWilsonFermion : public FermionAction<LatticeFermion,LatticeGaugeField>
class WilsonFermion5D : public FermionOperator<LatticeFermion,LatticeGaugeField>
{
public:
///////////////////////////////////////////////////////////////
@ -26,19 +26,21 @@ namespace Grid {
GridBase *FermionGrid(void) { return _FiveDimGrid;}
GridBase *FermionRedBlackGrid(void) { return _FiveDimRedBlackGrid;}
// override multiply
virtual RealD M (const LatticeFermion &in, LatticeFermion &out);
virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out);
// full checkerboard operations; leave unimplemented as abstract for now
//virtual RealD M (const LatticeFermion &in, LatticeFermion &out)=0;
//virtual RealD Mdag (const LatticeFermion &in, LatticeFermion &out)=0;
// 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);
// half checkerboard operations; leave unimplemented as abstract for now
// virtual void Meooe (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MeooeDag (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void Mooee (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeDag (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeInv (const LatticeFermion &in, LatticeFermion &out)=0;
// virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out)=0;
// non-hermitian hopping term; half cb or both
// 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);
@ -54,12 +56,12 @@ namespace Grid {
int dag);
// Constructors
FiveDimWilsonFermion(LatticeGaugeField &_Umu,
WilsonFermion5D(LatticeGaugeField &_Umu,
GridCartesian &FiveDimGrid,
GridRedBlackCartesian &FiveDimRedBlackGrid,
GridCartesian &FourDimGrid,
GridRedBlackCartesian &FourDimRedBlackGrid,
double _mass);
double _M5);
// DoubleStore
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@ -82,7 +84,6 @@ namespace Grid {
static const std::vector<int> displacements;
double M5;
double mass;
int Ls;
//Defines the stencils for even and odd

12
lib/stencil/Grid_stencil_common.cc
View File

@ -52,8 +52,8 @@ namespace Grid {
// up a table containing the npoint "neighbours" and whether they
// live in lattice or a comms buffer.
if ( !comm_dim ) {
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Local(point,dimension,shift,0x3);
@ -63,8 +63,8 @@ namespace Grid {
}
} else { // All permute extract done in comms phase prior to Stencil application
// So tables are the same whether comm_dim or splice_dim
sshift[0] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShift(_checkerboard,dimension,shift,Odd);
sshift[0] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Even);
sshift[1] = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,Odd);
if ( sshift[0] == sshift[1] ) {
Comms(point,dimension,shift,0x3);
} else {
@ -96,7 +96,7 @@ namespace Grid {
int cb= (cbmask==0x2)? Odd : Even;
int sshift = _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
int sshift = _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
int sx = (x+sshift)%rd;
int permute_slice=0;
@ -134,7 +134,7 @@ namespace Grid {
// send to one or more remote nodes.
int cb= (cbmask==0x2)? Odd : Even;
int sshift= _grid->CheckerBoardShift(_checkerboard,dimension,shift,cb);
int sshift= _grid->CheckerBoardShiftForCB(_checkerboard,dimension,shift,cb);
for(int x=0;x<rd;x++){