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

Better EO support letting Schur solver work

Browse Source
This commit is contained in:
Peter Boyle 2015-05-25 13:46:28 +01:00
parent 1d4b1c48cc
commit 201a110c51
3 changed files with 210 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
lib/qcd/Grid_qcd.h
View File

@ -10,9 +10,6 @@ namespace QCD {
static const int Nhs=2; // half spinor static const int Nhs=2; // half spinor
static const int Nds=8; // double stored gauge field static const int Nds=8; // double stored gauge field
static const int CbRed =0;
static const int CbBlack=1;
////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////
// QCD iMatrix types // QCD iMatrix types
// Index conventions: Lorentz x Spin x Colour // Index conventions: Lorentz x Spin x Colour

317
lib/qcd/Grid_qcd_wilson_dop.cc
View File

@ -4,8 +4,8 @@
namespace Grid { namespace Grid {
namespace QCD { namespace QCD {
const std::vector<int> WilsonMatrix::directions ({0,1,2,3, 0, 1, 2, 3,0}); const std::vector<int> WilsonMatrix::directions ({0,1,2,3, 0, 1, 2, 3});
const std::vector<int> WilsonMatrix::displacements({1,1,1,1,-1,-1,-1,-1,0}); const std::vector<int> WilsonMatrix::displacements({1,1,1,1,-1,-1,-1,-1});
// Should be in header? // Should be in header?
const int WilsonMatrix::Xp = 0; const int WilsonMatrix::Xp = 0;
@ -16,7 +16,6 @@ const int WilsonMatrix::Xm = 4;
const int WilsonMatrix::Ym = 5; const int WilsonMatrix::Ym = 5;
const int WilsonMatrix::Zm = 6; const int WilsonMatrix::Zm = 6;
const int WilsonMatrix::Tm = 7; const int WilsonMatrix::Tm = 7;
//const int WilsonMatrix::X0 = 8;
class WilsonCompressor { class WilsonCompressor {
public: public:
@ -72,20 +71,27 @@ const int WilsonMatrix::Tm = 7;
} }
}; };
WilsonMatrix::WilsonMatrix(LatticeGaugeField &_Umu,double _mass) WilsonMatrix::WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid, double _mass) :
: Stencil(Umu._grid,npoint,0,directions,displacements), CheckerBoardedSparseMatrixBase<LatticeFermion>(&Fgrid,&Hgrid),
Stencil ( _grid,npoint,Even,directions,displacements),
StencilEven(_cbgrid,npoint,Even,directions,displacements), // source is Even
StencilOdd (_cbgrid,npoint,Odd ,directions,displacements), // source is Odd
mass(_mass), mass(_mass),
Umu(_Umu._grid) Umu(_grid),
{ UmuEven(_cbgrid),
// Allocate the required comms buffer UmuOdd (_cbgrid)
grid = _Umu._grid; {
comm_buf.resize(Stencil._unified_buffer_size); // Allocate the required comms buffer
DoubleStore(Umu,_Umu); comm_buf.resize(Stencil._unified_buffer_size); // this is always big enough to contain EO
}
DoubleStore(Umu,_Umu);
pickCheckerboard(Even,UmuEven,Umu);
pickCheckerboard(Odd ,UmuOdd,Umu);
}
void WilsonMatrix::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu) void WilsonMatrix::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu)
{ {
LatticeColourMatrix U(grid); LatticeColourMatrix U(_grid);
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
U = peekIndex<LorentzIndex>(Umu,mu); U = peekIndex<LorentzIndex>(Umu,mu);
@ -97,47 +103,63 @@ void WilsonMatrix::DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeF
RealD WilsonMatrix::M(const LatticeFermion &in, LatticeFermion &out) RealD WilsonMatrix::M(const LatticeFermion &in, LatticeFermion &out)
{ {
this->Dhop(in,out,0); out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerNo);
out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun
return norm2(out); return norm2(out);
} }
RealD WilsonMatrix::Mdag(const LatticeFermion &in, LatticeFermion &out) RealD WilsonMatrix::Mdag(const LatticeFermion &in, LatticeFermion &out)
{ {
this->Dhop(in,out,1); out.checkerboard=in.checkerboard;
Dhop(in,out,DaggerYes);
out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun out = (4+mass)*in - 0.5*out ; // FIXME : axpby_norm! fusion fun
return norm2(out); return norm2(out);
} }
void WilsonMatrix::Meooe(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::Meooe(const LatticeFermion &in, LatticeFermion &out)
{ {
this->Dhop(in,out,0); if ( in.checkerboard == Odd ) {
out = 0.5*out; // FIXME : scale factor in Dhop DhopEO(in,out,DaggerNo);
} else {
DhopOE(in,out,DaggerNo);
}
out = (-0.5)*out; // FIXME : scale factor in Dhop
} }
void WilsonMatrix::MeooeDag(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::MeooeDag(const LatticeFermion &in, LatticeFermion &out)
{ {
this->Dhop(in,out,1); if ( in.checkerboard == Odd ) {
out = 0.5*out; // FIXME : scale factor in Dhop DhopEO(in,out,DaggerYes);
} else {
DhopOE(in,out,DaggerYes);
}
out = (-0.5)*out; // FIXME : scale factor in Dhop
} }
void WilsonMatrix::Mooee(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::Mooee(const LatticeFermion &in, LatticeFermion &out)
{ {
out.checkerboard = in.checkerboard;
out = (4.0+mass)*in; out = (4.0+mass)*in;
return ; return ;
} }
void WilsonMatrix::MooeeDag(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::MooeeDag(const LatticeFermion &in, LatticeFermion &out)
{ {
this->Mooee(in,out); out.checkerboard = in.checkerboard;
Mooee(in,out);
} }
void WilsonMatrix::MooeeInv(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::MooeeInv(const LatticeFermion &in, LatticeFermion &out)
{ {
out.checkerboard = in.checkerboard;
out = (1.0/(4.0+mass))*in; out = (1.0/(4.0+mass))*in;
return ; return ;
} }
void WilsonMatrix::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::MooeeInvDag(const LatticeFermion &in, LatticeFermion &out)
{ {
this->MooeeInv(in,out); out.checkerboard = in.checkerboard;
MooeeInv(in,out);
} }
void WilsonMatrix::DhopSite(int ss,const LatticeFermion &in, LatticeFermion &out) void WilsonMatrix::DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int ss,const LatticeFermion &in, LatticeFermion &out)
{ {
vHalfSpinColourVector tmp; vHalfSpinColourVector tmp;
vHalfSpinColourVector chi; vHalfSpinColourVector chi;
@ -145,79 +167,78 @@ void WilsonMatrix::DhopSite(int ss,const LatticeFermion &in, LatticeFermion &out
vHalfSpinColourVector Uchi; vHalfSpinColourVector Uchi;
int offset,local,perm, ptype; int offset,local,perm, ptype;
// int ss = Stencil._LebesgueReorder[sss]; //#define VERBOSE( A) if ( ss<10 ) { std::cout << "site " <<ss << " " #A " neigh " << offset << " perm "<< perm <<std::endl;}
int ssu= ss;
// Xp // Xp
offset = Stencil._offsets [Xp][ss]; offset = st._offsets [Xp][ss];
local = Stencil._is_local[Xp][ss]; local = st._is_local[Xp][ss];
perm = Stencil._permute[Xp][ss]; perm = st._permute[Xp][ss];
ptype = Stencil._permute_type[Xp]; ptype = st._permute_type[Xp];
if ( local && perm ) { if ( local && perm ) {
spProjXp(tmp,in._odata[offset]); spProjXp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjXp(chi,in._odata[offset]); spProjXp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Xp),&chi()); mult(&Uchi(),&U._odata[ss](Xp),&chi());
spReconXp(result,Uchi); spReconXp(result,Uchi);
// Yp // Yp
offset = Stencil._offsets [Yp][ss]; offset = st._offsets [Yp][ss];
local = Stencil._is_local[Yp][ss]; local = st._is_local[Yp][ss];
perm = Stencil._permute[Yp][ss]; perm = st._permute[Yp][ss];
ptype = Stencil._permute_type[Yp]; ptype = st._permute_type[Yp];
if ( local && perm ) { if ( local && perm ) {
spProjYp(tmp,in._odata[offset]); spProjYp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjYp(chi,in._odata[offset]); spProjYp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Yp),&chi()); mult(&Uchi(),&U._odata[ss](Yp),&chi());
accumReconYp(result,Uchi); accumReconYp(result,Uchi);
// Zp // Zp
offset = Stencil._offsets [Zp][ss]; offset = st._offsets [Zp][ss];
local = Stencil._is_local[Zp][ss]; local = st._is_local[Zp][ss];
perm = Stencil._permute[Zp][ss]; perm = st._permute[Zp][ss];
ptype = Stencil._permute_type[Zp]; ptype = st._permute_type[Zp];
if ( local && perm ) { if ( local && perm ) {
spProjZp(tmp,in._odata[offset]); spProjZp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjZp(chi,in._odata[offset]); spProjZp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Zp),&chi()); mult(&Uchi(),&U._odata[ss](Zp),&chi());
accumReconZp(result,Uchi); accumReconZp(result,Uchi);
// Tp // Tp
offset = Stencil._offsets [Tp][ss]; offset = st._offsets [Tp][ss];
local = Stencil._is_local[Tp][ss]; local = st._is_local[Tp][ss];
perm = Stencil._permute[Tp][ss]; perm = st._permute[Tp][ss];
ptype = Stencil._permute_type[Tp]; ptype = st._permute_type[Tp];
if ( local && perm ) { if ( local && perm ) {
spProjTp(tmp,in._odata[offset]); spProjTp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjTp(chi,in._odata[offset]); spProjTp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Tp),&chi()); mult(&Uchi(),&U._odata[ss](Tp),&chi());
accumReconTp(result,Uchi); accumReconTp(result,Uchi);
// Xm // Xm
offset = Stencil._offsets [Xm][ss]; offset = st._offsets [Xm][ss];
local = Stencil._is_local[Xm][ss]; local = st._is_local[Xm][ss];
perm = Stencil._permute[Xm][ss]; perm = st._permute[Xm][ss];
ptype = Stencil._permute_type[Xm]; ptype = st._permute_type[Xm];
if ( local && perm ) if ( local && perm )
{ {
@ -226,17 +247,17 @@ void WilsonMatrix::DhopSite(int ss,const LatticeFermion &in, LatticeFermion &out
} else if ( local ) { } else if ( local ) {
spProjXm(chi,in._odata[offset]); spProjXm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Xm),&chi()); mult(&Uchi(),&U._odata[ss](Xm),&chi());
accumReconXm(result,Uchi); accumReconXm(result,Uchi);
// Ym // Ym
offset = Stencil._offsets [Ym][ss]; offset = st._offsets [Ym][ss];
local = Stencil._is_local[Ym][ss]; local = st._is_local[Ym][ss];
perm = Stencil._permute[Ym][ss]; perm = st._permute[Ym][ss];
ptype = Stencil._permute_type[Ym]; ptype = st._permute_type[Ym];
if ( local && perm ) { if ( local && perm ) {
spProjYm(tmp,in._odata[offset]); spProjYm(tmp,in._odata[offset]);
@ -244,46 +265,49 @@ void WilsonMatrix::DhopSite(int ss,const LatticeFermion &in, LatticeFermion &out
} else if ( local ) { } else if ( local ) {
spProjYm(chi,in._odata[offset]); spProjYm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Ym),&chi()); mult(&Uchi(),&U._odata[ss](Ym),&chi());
accumReconYm(result,Uchi); accumReconYm(result,Uchi);
// Zm // Zm
offset = Stencil._offsets [Zm][ss]; offset = st._offsets [Zm][ss];
local = Stencil._is_local[Zm][ss]; local = st._is_local[Zm][ss];
perm = Stencil._permute[Zm][ss]; perm = st._permute[Zm][ss];
ptype = Stencil._permute_type[Zm]; ptype = st._permute_type[Zm];
if ( local && perm ) { if ( local && perm ) {
spProjZm(tmp,in._odata[offset]); spProjZm(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjZm(chi,in._odata[offset]); spProjZm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Zm),&chi()); mult(&Uchi(),&U._odata[ss](Zm),&chi());
accumReconZm(result,Uchi); accumReconZm(result,Uchi);
// Tm // Tm
offset = Stencil._offsets [Tm][ss]; offset = st._offsets [Tm][ss];
local = Stencil._is_local[Tm][ss]; local = st._is_local[Tm][ss];
perm = Stencil._permute[Tm][ss]; perm = st._permute[Tm][ss];
ptype = Stencil._permute_type[Tm]; ptype = st._permute_type[Tm];
if ( local && perm ) { if ( local && perm ) {
spProjTm(tmp,in._odata[offset]); spProjTm(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjTm(chi,in._odata[offset]); spProjTm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Tm),&chi()); mult(&Uchi(),&U._odata[ss](Tm),&chi());
accumReconTm(result,Uchi); accumReconTm(result,Uchi);
vstream(out._odata[ss],result); vstream(out._odata[ss],result);
} }
void WilsonMatrix::DhopSiteDag(int ss,const LatticeFermion &in, LatticeFermion &out)
void WilsonMatrix::DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int ss,const LatticeFermion &in, LatticeFermion &out)
{ {
vHalfSpinColourVector tmp; vHalfSpinColourVector tmp;
vHalfSpinColourVector chi; vHalfSpinColourVector chi;
@ -291,78 +315,76 @@ void WilsonMatrix::DhopSiteDag(int ss,const LatticeFermion &in, LatticeFermion &
vHalfSpinColourVector Uchi; vHalfSpinColourVector Uchi;
int offset,local,perm, ptype; int offset,local,perm, ptype;
int ssu= ss;
// Xp // Xp
offset = Stencil._offsets [Xm][ss]; offset = st._offsets [Xm][ss];
local = Stencil._is_local[Xm][ss]; local = st._is_local[Xm][ss];
perm = Stencil._permute[Xm][ss]; perm = st._permute[Xm][ss];
ptype = Stencil._permute_type[Xm]; ptype = st._permute_type[Xm];
if ( local && perm ) { if ( local && perm ) {
spProjXp(tmp,in._odata[offset]); spProjXp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjXp(chi,in._odata[offset]); spProjXp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Xm),&chi()); mult(&Uchi(),&U._odata[ss](Xm),&chi());
spReconXp(result,Uchi); spReconXp(result,Uchi);
// Yp // Yp
offset = Stencil._offsets [Ym][ss]; offset = st._offsets [Ym][ss];
local = Stencil._is_local[Ym][ss]; local = st._is_local[Ym][ss];
perm = Stencil._permute[Ym][ss]; perm = st._permute[Ym][ss];
ptype = Stencil._permute_type[Ym]; ptype = st._permute_type[Ym];
if ( local && perm ) { if ( local && perm ) {
spProjYp(tmp,in._odata[offset]); spProjYp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjYp(chi,in._odata[offset]); spProjYp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Ym),&chi()); mult(&Uchi(),&U._odata[ss](Ym),&chi());
accumReconYp(result,Uchi); accumReconYp(result,Uchi);
// Zp // Zp
offset = Stencil._offsets [Zm][ss]; offset = st._offsets [Zm][ss];
local = Stencil._is_local[Zm][ss]; local = st._is_local[Zm][ss];
perm = Stencil._permute[Zm][ss]; perm = st._permute[Zm][ss];
ptype = Stencil._permute_type[Zm]; ptype = st._permute_type[Zm];
if ( local && perm ) { if ( local && perm ) {
spProjZp(tmp,in._odata[offset]); spProjZp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjZp(chi,in._odata[offset]); spProjZp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Zm),&chi()); mult(&Uchi(),&U._odata[ss](Zm),&chi());
accumReconZp(result,Uchi); accumReconZp(result,Uchi);
// Tp // Tp
offset = Stencil._offsets [Tm][ss]; offset = st._offsets [Tm][ss];
local = Stencil._is_local[Tm][ss]; local = st._is_local[Tm][ss];
perm = Stencil._permute[Tm][ss]; perm = st._permute[Tm][ss];
ptype = Stencil._permute_type[Tm]; ptype = st._permute_type[Tm];
if ( local && perm ) { if ( local && perm ) {
spProjTp(tmp,in._odata[offset]); spProjTp(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjTp(chi,in._odata[offset]); spProjTp(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Tm),&chi()); mult(&Uchi(),&U._odata[ss](Tm),&chi());
accumReconTp(result,Uchi); accumReconTp(result,Uchi);
// Xm // Xm
offset = Stencil._offsets [Xp][ss]; offset = st._offsets [Xp][ss];
local = Stencil._is_local[Xp][ss]; local = st._is_local[Xp][ss];
perm = Stencil._permute[Xp][ss]; perm = st._permute[Xp][ss];
ptype = Stencil._permute_type[Xp]; ptype = st._permute_type[Xp];
if ( local && perm ) if ( local && perm )
{ {
@ -371,17 +393,16 @@ void WilsonMatrix::DhopSiteDag(int ss,const LatticeFermion &in, LatticeFermion &
} else if ( local ) { } else if ( local ) {
spProjXm(chi,in._odata[offset]); spProjXm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Xp),&chi()); mult(&Uchi(),&U._odata[ss](Xp),&chi());
accumReconXm(result,Uchi); accumReconXm(result,Uchi);
// Ym // Ym
offset = Stencil._offsets [Yp][ss]; offset = st._offsets [Yp][ss];
local = Stencil._is_local[Yp][ss]; local = st._is_local[Yp][ss];
perm = Stencil._permute[Yp][ss]; perm = st._permute[Yp][ss];
ptype = Stencil._permute_type[Yp]; ptype = st._permute_type[Yp];
if ( local && perm ) { if ( local && perm ) {
spProjYm(tmp,in._odata[offset]); spProjYm(tmp,in._odata[offset]);
@ -389,67 +410,97 @@ void WilsonMatrix::DhopSiteDag(int ss,const LatticeFermion &in, LatticeFermion &
} else if ( local ) { } else if ( local ) {
spProjYm(chi,in._odata[offset]); spProjYm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Yp),&chi()); mult(&Uchi(),&U._odata[ss](Yp),&chi());
accumReconYm(result,Uchi); accumReconYm(result,Uchi);
// Zm // Zm
offset = Stencil._offsets [Zp][ss]; offset = st._offsets [Zp][ss];
local = Stencil._is_local[Zp][ss]; local = st._is_local[Zp][ss];
perm = Stencil._permute[Zp][ss]; perm = st._permute[Zp][ss];
ptype = Stencil._permute_type[Zp]; ptype = st._permute_type[Zp];
if ( local && perm ) { if ( local && perm ) {
spProjZm(tmp,in._odata[offset]); spProjZm(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjZm(chi,in._odata[offset]); spProjZm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Zp),&chi()); mult(&Uchi(),&U._odata[ss](Zp),&chi());
accumReconZm(result,Uchi); accumReconZm(result,Uchi);
// Tm // Tm
offset = Stencil._offsets [Tp][ss]; offset = st._offsets [Tp][ss];
local = Stencil._is_local[Tp][ss]; local = st._is_local[Tp][ss];
perm = Stencil._permute[Tp][ss]; perm = st._permute[Tp][ss];
ptype = Stencil._permute_type[Tp]; ptype = st._permute_type[Tp];
if ( local && perm ) { if ( local && perm ) {
spProjTm(tmp,in._odata[offset]); spProjTm(tmp,in._odata[offset]);
permute(chi,tmp,ptype); permute(chi,tmp,ptype);
} else if ( local ) { } else if ( local ) {
spProjTm(chi,in._odata[offset]); spProjTm(chi,in._odata[offset]);
} else { } else {
chi=comm_buf[offset]; chi=buf[offset];
} }
mult(&Uchi(),&Umu._odata[ssu](Tp),&chi()); mult(&Uchi(),&U._odata[ss](Tp),&chi());
accumReconTm(result,Uchi); accumReconTm(result,Uchi);
vstream(out._odata[ss],result); vstream(out._odata[ss],result);
} }
void WilsonMatrix::Dhop(const LatticeFermion &in, LatticeFermion &out,int dag) void WilsonMatrix::DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField & U,
const LatticeFermion &in, LatticeFermion &out,int dag)
{ {
assert((dag==0) ||(dag==1)); assert((dag==DaggerNo) ||(dag==DaggerYes));
WilsonCompressor compressor(dag); WilsonCompressor compressor(dag);
Stencil.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
if ( dag ) { st.HaloExchange<vSpinColourVector,vHalfSpinColourVector,WilsonCompressor>(in,comm_buf,compressor);
if ( dag == DaggerYes ) {
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<grid->oSites();sss++){ for(int sss=0;sss<in._grid->oSites();sss++){
DhopSiteDag(sss,in,out); DhopSiteDag(st,U,comm_buf,sss,in,out);
} }
} else { } else {
PARALLEL_FOR_LOOP PARALLEL_FOR_LOOP
for(int sss=0;sss<grid->oSites();sss++){ for(int sss=0;sss<in._grid->oSites();sss++){
DhopSite(sss,in,out); DhopSite(st,U,comm_buf,sss,in,out);
} }
} }
} }
void WilsonMatrix::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 WilsonMatrix::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 WilsonMatrix::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);
} }
}}

36
lib/qcd/Grid_qcd_wilson_dop.h
View File

@ -11,14 +11,20 @@ namespace Grid {
//NB r=1; //NB r=1;
public: public:
double mass; double mass;
GridBase *grid; // GridBase * grid; // Inherited
// GridBase * cbgrid;
// Copy of the gauge field //Defines the stencils for even and odd
LatticeDoubledGaugeField Umu; CartesianStencil Stencil;
CartesianStencil StencilEven;
CartesianStencil StencilOdd;
//Defines the stencil // Copy of the gauge field , with even and odd subsets
CartesianStencil Stencil; LatticeDoubledGaugeField Umu;
static const int npoint=9; LatticeDoubledGaugeField UmuEven;
LatticeDoubledGaugeField UmuOdd;
static const int npoint=8;
static const std::vector<int> directions ; static const std::vector<int> directions ;
static const std::vector<int> displacements; static const std::vector<int> displacements;
static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm; static const int Xp,Xm,Yp,Ym,Zp,Zm,Tp,Tm;
@ -27,7 +33,7 @@ namespace Grid {
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > comm_buf; std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > comm_buf;
// Constructor // Constructor
WilsonMatrix(LatticeGaugeField &Umu,double mass); WilsonMatrix(LatticeGaugeField &_Umu,GridCartesian &Fgrid,GridRedBlackCartesian &Hgrid,double _mass);
// DoubleStore // DoubleStore
void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu); void DoubleStore(LatticeDoubledGaugeField &Uds,const LatticeGaugeField &Umu);
@ -45,9 +51,19 @@ namespace Grid {
virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out); virtual void MooeeInvDag (const LatticeFermion &in, LatticeFermion &out);
// non-hermitian hopping term; half cb or both // non-hermitian hopping term; half cb or both
void Dhop(const LatticeFermion &in, LatticeFermion &out,int dag); void Dhop (const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopSite (int ss,const LatticeFermion &in, LatticeFermion &out); void DhopOE(const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopSiteDag(int ss,const LatticeFermion &in, LatticeFermion &out); void DhopEO(const LatticeFermion &in, LatticeFermion &out,int dag);
void DhopInternal(CartesianStencil & st,LatticeDoubledGaugeField &U,
const LatticeFermion &in, LatticeFermion &out,int dag);
// These ones will need to be package intelligently. WilsonType base class
// for use by DWF etc..
void DhopSite(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int ss,const LatticeFermion &in, LatticeFermion &out);
void DhopSiteDag(CartesianStencil &st,LatticeDoubledGaugeField &U,
std::vector<vHalfSpinColourVector,alignedAllocator<vHalfSpinColourVector> > &buf,
int ss,const LatticeFermion &in, LatticeFermion &out);
typedef iScalar<iMatrix<vComplex, Nc> > matrix; typedef iScalar<iMatrix<vComplex, Nc> > matrix;