portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-04-09 21:50:45 +01:00
Code Releases Activity

Zero changes, acceleartor on kernels and some thread loop changes

Browse Source
This commit is contained in:
paboyle 2018-01-27 23:47:38 +00:00
parent 45df59720e
commit 2d0bcc2606
40 changed files with 174 additions and 175 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
lib/qcd/action/fermion/CayleyFermion5Ddense.cc
View File

@ -96,15 +96,14 @@ void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField
} }
// For the non-vectorised s-direction this is simple // For the non-vectorised s-direction this is simple
thread_loop( (auto site=0;site<vol;site++), {
for(auto site=0;site<vol;site++){
SiteSpinor SiteChi; SiteSpinor SiteChi;
SiteHalfSpinor SitePplus; SiteHalfSpinor SitePplus;
SiteHalfSpinor SitePminus; SiteHalfSpinor SitePminus;
for(int s1=0;s1<Ls;s1++){ for(int s1=0;s1<Ls;s1++){
SiteChi =zero; SiteChi =Zero();
for(int s2=0;s2<Ls;s2++){ for(int s2=0;s2<Ls;s2++){
int lex2 = s2+Ls*site; int lex2 = s2+Ls*site;
@ -120,7 +119,7 @@ void CayleyFermion5D<Impl>::MooeeInternal(const FermionField &psi, FermionField
} }
chi[s1+Ls*site] = SiteChi*0.5; chi[s1+Ls*site] = SiteChi*0.5;
} }
} });
} }
#ifdef CAYLEY_DPERP_DENSE #ifdef CAYLEY_DPERP_DENSE

8
lib/qcd/action/fermion/CayleyFermion5Dvec.cc
View File

@ -360,8 +360,8 @@ void CayleyFermion5D<Impl>::MooeeInternalAsm(const FermionField &psi, FermionFie
int lex=s2+LLs*site; int lex=s2+LLs*site;
if ( s2==0 && l==0) { if ( s2==0 && l==0) {
SiteChiP=zero; SiteChiP=Zero();
SiteChiM=zero; SiteChiM=Zero();
} }
for(int sp=0;sp<2;sp++){ for(int sp=0;sp<2;sp++){
@ -532,8 +532,8 @@ void CayleyFermion5D<Impl>::MooeeInternalZAsm(const FermionField &psi, FermionFi
int lex=s2+LLs*site; int lex=s2+LLs*site;
if ( s2==0 && l==0) { if ( s2==0 && l==0) {
SiteChiP=zero; SiteChiP=Zero();
SiteChiM=zero; SiteChiM=Zero();
} }
for(int sp=0;sp<2;sp++){ for(int sp=0;sp<2;sp++){

2
lib/qcd/action/fermion/DomainWallEOFAFermion.cc
View File

@ -69,7 +69,7 @@ void DomainWallEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& D
{ {
int Ls = this->Ls; int Ls = this->Ls;
Din = zero; Din = Zero();
if((sign == 1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, Ls-1, 0); } if((sign == 1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, Ls-1, 0); }
else if((sign == -1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); } else if((sign == -1) && (dag == 0)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, 0); }
else if((sign == 1 ) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, Ls-1); } else if((sign == 1 ) && (dag == 1)){ axpby_ssp(Din, 0.0, psi, 1.0, psi, 0, Ls-1); }

2
lib/qcd/action/fermion/DomainWallEOFAFermiondense.cc
View File

@ -106,7 +106,7 @@ void DomainWallEOFAFermion<Impl>::MooeeInternal(const FermionField& psi, Fermion
SiteHalfSpinor SitePminus; SiteHalfSpinor SitePminus;
for(int s1=0; s1<Ls; s1++){ for(int s1=0; s1<Ls; s1++){
SiteChi = zero; SiteChi = Zero();
for(int s2=0; s2<Ls; s2++){ for(int s2=0; s2<Ls; s2++){
int lex2 = s2 + Ls*site; int lex2 = s2 + Ls*site;
if(PplusMat(s1,s2) != 0.0){ if(PplusMat(s1,s2) != 0.0){

4
lib/qcd/action/fermion/DomainWallEOFAFermionvec.cc
View File

@ -362,8 +362,8 @@ void DomainWallEOFAFermion<Impl>::MooeeInternalAsm(const FermionField& psi, Ferm
int lex = s2 + LLs*site; int lex = s2 + LLs*site;
if( s2==0 && l==0 ){ if( s2==0 && l==0 ){
SiteChiP=zero; SiteChiP=Zero();
SiteChiM=zero; SiteChiM=Zero();
} }
for(int sp=0; sp<2; sp++){ for(int sp=0; sp<2; sp++){

4
lib/qcd/action/fermion/FermionOperator.h
View File

@ -81,8 +81,8 @@ public:
virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);}; virtual void MDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDeriv(mat,U,V,dag);};
virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);}; virtual void MoeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivOE(mat,U,V,dag);};
virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);}; virtual void MeoDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){DhopDerivEO(mat,U,V,dag);};
virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; // Clover can override these virtual void MooDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=Zero();}; // Clover can override these
virtual void MeeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=zero;}; virtual void MeeDeriv(GaugeField &mat,const FermionField &U,const FermionField &V,int dag){mat=Zero();};
virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDeriv (GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;
virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0; virtual void DhopDerivEO(GaugeField &mat,const FermionField &U,const FermionField &V,int dag)=0;

8
lib/qcd/action/fermion/FermionOperatorImpl.h
View File

@ -266,7 +266,7 @@ public:
int Ls=Btilde.Grid()->_fdimensions[0]; int Ls=Btilde.Grid()->_fdimensions[0];
GaugeLinkField tmp(mat.Grid()); GaugeLinkField tmp(mat.Grid());
tmp = zero; tmp = Zero();
parallel_for(int sss=0;sss<tmp.Grid()->oSites();sss++){ parallel_for(int sss=0;sss<tmp.Grid()->oSites();sss++){
int sU=sss; int sU=sss;
@ -406,7 +406,7 @@ public:
unsigned int dimU = grid->Nd(); unsigned int dimU = grid->Nd();
unsigned int dimF = Bgrid->Nd(); unsigned int dimF = Bgrid->Nd();
GaugeLinkField tmp(grid); GaugeLinkField tmp(grid);
tmp = zero; tmp = Zero();
// FIXME // FIXME
// Current implementation works, thread safe, probably suboptimal // Current implementation works, thread safe, probably suboptimal
@ -417,7 +417,7 @@ public:
std::vector<typename result_type::scalar_object> vres(Bgrid->Nsimd()); std::vector<typename result_type::scalar_object> vres(Bgrid->Nsimd());
std::vector<int> ocoor; grid->oCoorFromOindex(ocoor,so); std::vector<int> ocoor; grid->oCoorFromOindex(ocoor,so);
for (int si = 0; si < tmp.Grid()->iSites(); si++){ for (int si = 0; si < tmp.Grid()->iSites(); si++){
typename result_type::scalar_object scalar_object; scalar_object = zero; typename result_type::scalar_object scalar_object; scalar_object = Zero();
std::vector<int> local_coor; std::vector<int> local_coor;
std::vector<int> icoor; grid->iCoorFromIindex(icoor,si); std::vector<int> icoor; grid->iCoorFromIindex(icoor,si);
grid->InOutCoorToLocalCoor(ocoor, icoor, local_coor); grid->InOutCoorToLocalCoor(ocoor, icoor, local_coor);
@ -639,7 +639,7 @@ public:
int Ls = Btilde.Grid()->_fdimensions[0]; int Ls = Btilde.Grid()->_fdimensions[0];
GaugeLinkField tmp(mat.Grid()); GaugeLinkField tmp(mat.Grid());
tmp = zero; tmp = Zero();
parallel_for(int ss = 0; ss < tmp.Grid()->oSites(); ss++) { parallel_for(int ss = 0; ss < tmp.Grid()->oSites(); ss++) {
for (int s = 0; s < Ls; s++) { for (int s = 0; s < Ls; s++) {
int sF = s + Ls * ss; int sF = s + Ls * ss;

2
lib/qcd/action/fermion/MobiusEOFAFermion.cc
View File

@ -87,7 +87,7 @@ void MobiusEOFAFermion<Impl>::Omega(const FermionField& psi, FermionField& Din,
int Ls = this->Ls; int Ls = this->Ls;
RealD alpha = this->alpha; RealD alpha = this->alpha;
Din = zero; Din = Zero();
if((sign == 1) && (dag == 0)) { // \Omega_{+} if((sign == 1) && (dag == 0)) { // \Omega_{+}
for(int s=0; s<Ls; ++s){ for(int s=0; s<Ls; ++s){
axpby_ssp(Din, 0.0, psi, 2.0*std::pow(1.0-alpha,Ls-s-1)/std::pow(1.0+alpha,Ls-s), psi, s, 0); axpby_ssp(Din, 0.0, psi, 2.0*std::pow(1.0-alpha,Ls-s-1)/std::pow(1.0+alpha,Ls-s), psi, s, 0);

2
lib/qcd/action/fermion/MobiusEOFAFermioncache.cc
View File

@ -175,7 +175,7 @@ void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi, const Fermio
this->M5Dtime -= usecond(); this->M5Dtime -= usecond();
parallel_for(int ss=0; ss<grid->oSites(); ss+=Ls){ parallel_for(int ss=0; ss<grid->oSites(); ss+=Ls){
chi[ss+Ls-1] = zero; chi[ss+Ls-1] = Zero();
auto tmp = psi[0]; auto tmp = psi[0];
for(int s=0; s<Ls; s++){ for(int s=0; s<Ls; s++){
if(s==0) { if(s==0) {

2
lib/qcd/action/fermion/MobiusEOFAFermiondense.cc
View File

@ -131,7 +131,7 @@ void MobiusEOFAFermion<Impl>::MooeeInternal(const FermionField& psi, FermionFiel
SiteHalfSpinor SitePminus; SiteHalfSpinor SitePminus;
for(int s1=0; s1<Ls; s1++){ for(int s1=0; s1<Ls; s1++){
SiteChi = zero; SiteChi = Zero();
for(int s2=0; s2<Ls; s2++){ for(int s2=0; s2<Ls; s2++){
int lex2 = s2 + Ls*site; int lex2 = s2 + Ls*site;
if(PplusMat(s1,s2) != 0.0){ if(PplusMat(s1,s2) != 0.0){

4
lib/qcd/action/fermion/MobiusEOFAFermionvec.cc
View File

@ -737,8 +737,8 @@ void MobiusEOFAFermion<Impl>::MooeeInternalAsm(const FermionField& psi, FermionF
int lex = s2 + LLs*site; int lex = s2 + LLs*site;
if( s2==0 && l==0 ){ if( s2==0 && l==0 ){
SiteChiP=zero; SiteChiP=Zero();
SiteChiM=zero; SiteChiM=Zero();
} }
for(int sp=0; sp<2; sp++){ for(int sp=0; sp<2; sp++){

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

@ -152,11 +152,11 @@ void WilsonFermion<Impl>::MomentumSpacePropagator(FermionField &out, const Fermi
std::vector<int> latt_size = _grid->_fdimensions; std::vector<int> latt_size = _grid->_fdimensions;
FermionField num (_grid); num = zero; FermionField num (_grid); num = Zero();
LatComplex wilson(_grid); wilson= zero; LatComplex wilson(_grid); wilson= Zero();
LatComplex one (_grid); one = ScalComplex(1.0,0.0); LatComplex one (_grid); one = ScalComplex(1.0,0.0);
LatComplex denom(_grid); denom= zero; LatComplex denom(_grid); denom= Zero();
LatComplex kmu(_grid); LatComplex kmu(_grid);
ScalComplex ci(0.0,1.0); ScalComplex ci(0.0,1.0);
// momphase = n * 2pi / L // momphase = n * 2pi / L
@ -360,7 +360,7 @@ void WilsonFermion<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
conformable(_grid, q_in_2.Grid()); conformable(_grid, q_in_2.Grid());
conformable(_grid, q_out.Grid()); conformable(_grid, q_out.Grid());
PropagatorField tmp1(_grid), tmp2(_grid); PropagatorField tmp1(_grid), tmp2(_grid);
q_out = zero; q_out = Zero();
// Forward, need q1(x + mu), q2(x). Backward, need q1(x), q2(x + mu). // Forward, need q1(x + mu), q2(x). Backward, need q1(x), q2(x + mu).
// Inefficient comms method but not performance critical. // Inefficient comms method but not performance critical.
@ -397,7 +397,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
unsigned int LLt = GridDefaultLatt()[Tp]; unsigned int LLt = GridDefaultLatt()[Tp];
// Momentum projection // Momentum projection
ph = zero; ph = Zero();
for(unsigned int mu = 0; mu < Nd - 1; mu++) for(unsigned int mu = 0; mu < Nd - 1; mu++)
{ {
LatticeCoordinate(coor, mu); LatticeCoordinate(coor, mu);
@ -405,7 +405,7 @@ void WilsonFermion<Impl>::SeqConservedCurrent(PropagatorField &q_in,
} }
ph = exp((Real)(2*M_PI)*i*ph); ph = exp((Real)(2*M_PI)*i*ph);
q_out = zero; q_out = Zero();
LatticeInteger coords(_grid); LatticeInteger coords(_grid);
LatticeCoordinate(coords, Tp); LatticeCoordinate(coords, Tp);

30
lib/qcd/action/fermion/WilsonFermion5D.cc
View File

@ -583,14 +583,14 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHt(FermionField &out,const Fe
std::vector<int> latt_size = _grid->_fdimensions; std::vector<int> latt_size = _grid->_fdimensions;
FermionField num (_grid); num = zero; FermionField num (_grid); num = Zero();
LatComplex sk(_grid); sk = zero; LatComplex sk(_grid); sk = Zero();
LatComplex sk2(_grid); sk2= zero; LatComplex sk2(_grid); sk2= Zero();
LatComplex W(_grid); W= zero; LatComplex W(_grid); W= Zero();
LatComplex a(_grid); a= zero; LatComplex a(_grid); a= Zero();
LatComplex one (_grid); one = ScalComplex(1.0,0.0); LatComplex one (_grid); one = ScalComplex(1.0,0.0);
LatComplex denom(_grid); denom= zero; LatComplex denom(_grid); denom= Zero();
LatComplex cosha(_grid); LatComplex cosha(_grid);
LatComplex kmu(_grid); LatComplex kmu(_grid);
LatComplex Wea(_grid); LatComplex Wea(_grid);
@ -661,16 +661,16 @@ void WilsonFermion5D<Impl>::MomentumSpacePropagatorHw(FermionField &out,const Fe
std::vector<int> latt_size = _grid->_fdimensions; std::vector<int> latt_size = _grid->_fdimensions;
LatComplex sk(_grid); sk = zero; LatComplex sk(_grid); sk = Zero();
LatComplex sk2(_grid); sk2= zero; LatComplex sk2(_grid); sk2= Zero();
LatComplex w_k(_grid); w_k= zero; LatComplex w_k(_grid); w_k= Zero();
LatComplex b_k(_grid); b_k= zero; LatComplex b_k(_grid); b_k= Zero();
LatComplex one (_grid); one = ScalComplex(1.0,0.0); LatComplex one (_grid); one = ScalComplex(1.0,0.0);
FermionField num (_grid); num = zero; FermionField num (_grid); num = Zero();
LatComplex denom(_grid); denom= zero; LatComplex denom(_grid); denom= Zero();
LatComplex kmu(_grid); LatComplex kmu(_grid);
ScalComplex ci(0.0,1.0); ScalComplex ci(0.0,1.0);
@ -733,7 +733,7 @@ void WilsonFermion5D<Impl>::ContractConservedCurrent(PropagatorField &q_in_1,
conformable(_FourDimGrid, q_out.Grid()); conformable(_FourDimGrid, q_out.Grid());
PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid()); PropagatorField tmp1(FermionGrid()), tmp2(FermionGrid());
unsigned int LLs = q_in_1.Grid()->_rdimensions[0]; unsigned int LLs = q_in_1.Grid()->_rdimensions[0];
q_out = zero; q_out = Zero();
// Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s), // Forward, need q1(x + mu, s), q2(x, Ls - 1 - s). Backward, need q1(x, s),
// q2(x + mu, Ls - 1 - s). 5D lattice so shift 4D coordinate mu by one. // q2(x + mu, Ls - 1 - s). 5D lattice so shift 4D coordinate mu by one.
@ -797,7 +797,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
unsigned int LLt = GridDefaultLatt()[Tp]; unsigned int LLt = GridDefaultLatt()[Tp];
// Momentum projection. // Momentum projection.
ph = zero; ph = Zero();
for(unsigned int nu = 0; nu < Nd - 1; nu++) for(unsigned int nu = 0; nu < Nd - 1; nu++)
{ {
// Shift coordinate lattice index by 1 to account for 5th dimension. // Shift coordinate lattice index by 1 to account for 5th dimension.
@ -806,7 +806,7 @@ void WilsonFermion5D<Impl>::SeqConservedCurrent(PropagatorField &q_in,
} }
ph = exp((Real)(2*M_PI)*i*ph); ph = exp((Real)(2*M_PI)*i*ph);
q_out = zero; q_out = Zero();
LatticeInteger coords(_FourDimGrid); LatticeInteger coords(_FourDimGrid);
LatticeCoordinate(coords, Tp); LatticeCoordinate(coords, Tp);

38
lib/qcd/action/fermion/WilsonKernels.cc
View File

@ -36,7 +36,7 @@ int WilsonKernelsStatic::Opt = WilsonKernelsStatic::OptGeneric;
int WilsonKernelsStatic::Comms = WilsonKernelsStatic::CommsAndCompute; int WilsonKernelsStatic::Comms = WilsonKernelsStatic::CommsAndCompute;
template <class Impl> template <class Impl>
WilsonKernels<Impl>::WilsonKernels(const ImplParams &p) : Base(p){}; accelerator WilsonKernels<Impl>::WilsonKernels(const ImplParams &p) : Base(p){};
//////////////////////////////////////////// ////////////////////////////////////////////
// Generic implementation; move to different file? // Generic implementation; move to different file?
@ -103,9 +103,9 @@ WilsonKernels<Impl>::WilsonKernels(const ImplParams &p) : Base(p){};
// All legs kernels ; comms then compute // All legs kernels ; comms then compute
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSiteDag(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
SiteHalfSpinor tmp; SiteHalfSpinor tmp;
SiteHalfSpinor chi; SiteHalfSpinor chi;
@ -127,9 +127,9 @@ void WilsonKernels<Impl>::GenericDhopSiteDag(StencilImpl &st, LebesgueOrder &lo,
}; };
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
SiteHalfSpinor tmp; SiteHalfSpinor tmp;
SiteHalfSpinor chi; SiteHalfSpinor chi;
@ -153,7 +153,7 @@ void WilsonKernels<Impl>::GenericDhopSite(StencilImpl &st, LebesgueOrder &lo, Do
// Interior kernels // Interior kernels
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
@ -165,7 +165,7 @@ void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilImpl &st, LebesgueOrder &
StencilEntry *SE; StencilEntry *SE;
int ptype; int ptype;
result=zero; result=Zero();
GENERIC_STENCIL_LEG_INT(Xp,spProjXp,accumReconXp); GENERIC_STENCIL_LEG_INT(Xp,spProjXp,accumReconXp);
GENERIC_STENCIL_LEG_INT(Yp,spProjYp,accumReconYp); GENERIC_STENCIL_LEG_INT(Yp,spProjYp,accumReconYp);
GENERIC_STENCIL_LEG_INT(Zp,spProjZp,accumReconZp); GENERIC_STENCIL_LEG_INT(Zp,spProjZp,accumReconZp);
@ -178,9 +178,9 @@ void WilsonKernels<Impl>::GenericDhopSiteDagInt(StencilImpl &st, LebesgueOrder &
}; };
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSiteInt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSiteInt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
SiteHalfSpinor tmp; SiteHalfSpinor tmp;
SiteHalfSpinor chi; SiteHalfSpinor chi;
@ -189,7 +189,7 @@ void WilsonKernels<Impl>::GenericDhopSiteInt(StencilImpl &st, LebesgueOrder &lo,
SiteSpinor result; SiteSpinor result;
StencilEntry *SE; StencilEntry *SE;
int ptype; int ptype;
result=zero; result=Zero();
GENERIC_STENCIL_LEG_INT(Xm,spProjXp,accumReconXp); GENERIC_STENCIL_LEG_INT(Xm,spProjXp,accumReconXp);
GENERIC_STENCIL_LEG_INT(Ym,spProjYp,accumReconYp); GENERIC_STENCIL_LEG_INT(Ym,spProjYp,accumReconYp);
GENERIC_STENCIL_LEG_INT(Zm,spProjZp,accumReconZp); GENERIC_STENCIL_LEG_INT(Zm,spProjZp,accumReconZp);
@ -204,7 +204,7 @@ void WilsonKernels<Impl>::GenericDhopSiteInt(StencilImpl &st, LebesgueOrder &lo,
// Exterior kernels // Exterior kernels
//////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
@ -216,7 +216,7 @@ void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilImpl &st, LebesgueOrder &
StencilEntry *SE; StencilEntry *SE;
int ptype; int ptype;
int nmu=0; int nmu=0;
result=zero; result=Zero();
GENERIC_STENCIL_LEG_EXT(Xp,spProjXp,accumReconXp); GENERIC_STENCIL_LEG_EXT(Xp,spProjXp,accumReconXp);
GENERIC_STENCIL_LEG_EXT(Yp,spProjYp,accumReconYp); GENERIC_STENCIL_LEG_EXT(Yp,spProjYp,accumReconYp);
GENERIC_STENCIL_LEG_EXT(Zp,spProjZp,accumReconZp); GENERIC_STENCIL_LEG_EXT(Zp,spProjZp,accumReconZp);
@ -231,7 +231,7 @@ void WilsonKernels<Impl>::GenericDhopSiteDagExt(StencilImpl &st, LebesgueOrder &
}; };
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::GenericDhopSiteExt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U, accelerator void WilsonKernels<Impl>::GenericDhopSiteExt(StencilImpl &st, LebesgueOrder &lo, DoubledGaugeField &U,
SiteHalfSpinor *buf, int sF, SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out) int sU, const FermionField &in, FermionField &out)
{ {
@ -243,7 +243,7 @@ void WilsonKernels<Impl>::GenericDhopSiteExt(StencilImpl &st, LebesgueOrder &lo,
StencilEntry *SE; StencilEntry *SE;
int ptype; int ptype;
int nmu=0; int nmu=0;
result=zero; result=Zero();
GENERIC_STENCIL_LEG_EXT(Xm,spProjXp,accumReconXp); GENERIC_STENCIL_LEG_EXT(Xm,spProjXp,accumReconXp);
GENERIC_STENCIL_LEG_EXT(Ym,spProjYp,accumReconYp); GENERIC_STENCIL_LEG_EXT(Ym,spProjYp,accumReconYp);
GENERIC_STENCIL_LEG_EXT(Zm,spProjZp,accumReconZp); GENERIC_STENCIL_LEG_EXT(Zm,spProjZp,accumReconZp);
@ -258,8 +258,8 @@ void WilsonKernels<Impl>::GenericDhopSiteExt(StencilImpl &st, LebesgueOrder &lo,
}; };
template <class Impl> template <class Impl>
void WilsonKernels<Impl>::DhopDirK( StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor *buf, int sF, accelerator void WilsonKernels<Impl>::DhopDirK( StencilImpl &st, DoubledGaugeField &U,SiteHalfSpinor *buf, int sF,
int sU, const FermionField &in, FermionField &out, int dir, int gamma) { int sU, const FermionField &in, FermionField &out, int dir, int gamma) {
SiteHalfSpinor tmp; SiteHalfSpinor tmp;
SiteHalfSpinor chi; SiteHalfSpinor chi;

24
lib/qcd/action/fermion/WilsonKernelsHand.cc
View File

@ -544,18 +544,18 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
Simd U_21; Simd U_21;
#define ZERO_RESULT \ #define ZERO_RESULT \
result_00=zero; \ result_00=Zero(); \
result_01=zero; \ result_01=Zero(); \
result_02=zero; \ result_02=Zero(); \
result_10=zero; \ result_10=Zero(); \
result_11=zero; \ result_11=Zero(); \
result_12=zero; \ result_12=Zero(); \
result_20=zero; \ result_20=Zero(); \
result_21=zero; \ result_21=Zero(); \
result_22=zero; \ result_22=Zero(); \
result_30=zero; \ result_30=Zero(); \
result_31=zero; \ result_31=Zero(); \
result_32=zero; result_32=Zero();
#define Chimu_00 Chi_00 #define Chimu_00 Chi_00
#define Chimu_01 Chi_01 #define Chimu_01 Chi_01

4
lib/qcd/action/gauge/GaugeImplTypes.h
View File

@ -88,7 +88,7 @@ public:
static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) { static inline void generate_momenta(Field &P, GridParallelRNG &pRNG) {
// specific for SU gauge fields // specific for SU gauge fields
LinkField Pmu(P.Grid()); LinkField Pmu(P.Grid());
Pmu = zero; Pmu = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); SU<Nrepresentation>::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu);
PokeIndex<LorentzIndex>(P, Pmu, mu); PokeIndex<LorentzIndex>(P, Pmu, mu);
@ -113,7 +113,7 @@ public:
static inline RealD FieldSquareNorm(Field& U){ static inline RealD FieldSquareNorm(Field& U){
LatticeComplex Hloc(U.Grid()); LatticeComplex Hloc(U.Grid());
Hloc = zero; Hloc = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
auto Umu = PeekIndex<LorentzIndex>(U, mu); auto Umu = PeekIndex<LorentzIndex>(U, mu);
Hloc += trace(Umu * Umu); Hloc += trace(Umu * Umu);

8
lib/qcd/action/gauge/Photon.h
View File

@ -108,7 +108,7 @@ void Photon<Gimpl>::invKHatSquared(GaugeLinkField &out)
TComplex Tzero= Complex(0.0,0.0); TComplex Tzero= Complex(0.0,0.0);
one = Complex(1.0,0.0); one = Complex(1.0,0.0);
out = zero; out = Zero();
for(int mu = 0; mu < nd; mu++) for(int mu = 0; mu < nd; mu++)
{ {
Real twoPiL = M_PI*2./l[mu]; Real twoPiL = M_PI*2./l[mu];
@ -144,7 +144,7 @@ void Photon<Gimpl>::zmSub(GaugeLinkField &out)
LatticeInteger spNrm(grid), coor(grid); LatticeInteger spNrm(grid), coor(grid);
GaugeLinkField z(grid); GaugeLinkField z(grid);
spNrm = zero; spNrm = Zero();
for(int d = 0; d < grid->_ndimension - 1; d++) for(int d = 0; d < grid->_ndimension - 1; d++)
{ {
LatticeCoordinate(coor,d); LatticeCoordinate(coor,d);
@ -252,7 +252,7 @@ void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng,
// //
// std::vector<int> latt_size = grid->_fdimensions; // std::vector<int> latt_size = grid->_fdimensions;
// //
// LatComplex denom(grid); denom= zero; // LatComplex denom(grid); denom= Zero();
// LatComplex one(grid); one = ScalComplex(1.0,0.0); // LatComplex one(grid); one = ScalComplex(1.0,0.0);
// LatComplex kmu(grid); // LatComplex kmu(grid);
// //
@ -278,7 +278,7 @@ void Photon<Gimpl>::StochasticField(GaugeField &out, GridParallelRNG &rng,
// //
// pokeSite(Tzero,denom,zero_mode); // pokeSite(Tzero,denom,zero_mode);
// //
// out = zero; // out = Zero();
// out = in*denom; // out = in*denom;
// }; // };

16
lib/qcd/action/pseudofermion/ExactOneFlavourRatio.h
View File

@ -131,7 +131,7 @@ public:
spProj(eta, tmp[0], -1, Lop.Ls); spProj(eta, tmp[0], -1, Lop.Ls);
Lop.Omega(tmp[0], tmp[1], -1, 0); Lop.Omega(tmp[0], tmp[1], -1, 0);
G5R5(CG_src, tmp[1]); G5R5(CG_src, tmp[1]);
tmp[1] = zero; tmp[1] = Zero();
for(int k=0; k<param.degree; ++k){ for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] ); gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
Lop.RefreshShiftCoefficients(-gamma_l); Lop.RefreshShiftCoefficients(-gamma_l);
@ -141,7 +141,7 @@ public:
Solver(Lop, CG_src, CG_soln); Solver(Lop, CG_src, CG_soln);
prev_solns.push_back(CG_soln); prev_solns.push_back(CG_soln);
} else { } else {
CG_soln = zero; // Just use zero as the initial guess CG_soln = Zero(); // Just use zero as the initial guess
Solver(Lop, CG_src, CG_soln); Solver(Lop, CG_src, CG_soln);
} }
Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back Lop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
@ -157,7 +157,7 @@ public:
spProj(eta, tmp[0], 1, Rop.Ls); spProj(eta, tmp[0], 1, Rop.Ls);
Rop.Omega(tmp[0], tmp[1], 1, 0); Rop.Omega(tmp[0], tmp[1], 1, 0);
G5R5(CG_src, tmp[1]); G5R5(CG_src, tmp[1]);
tmp[1] = zero; tmp[1] = Zero();
if(use_heatbath_forecasting){ prev_solns.clear(); } // empirically, LH solns don't help for RH solves if(use_heatbath_forecasting){ prev_solns.clear(); } // empirically, LH solns don't help for RH solves
for(int k=0; k<param.degree; ++k){ for(int k=0; k<param.degree; ++k){
gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] ); gamma_l = 1.0 / ( 1.0 + PowerNegHalf.poles[k] );
@ -168,7 +168,7 @@ public:
Solver(Rop, CG_src, CG_soln); Solver(Rop, CG_src, CG_soln);
prev_solns.push_back(CG_soln); prev_solns.push_back(CG_soln);
} else { } else {
CG_soln = zero; CG_soln = Zero();
Solver(Rop, CG_src, CG_soln); Solver(Rop, CG_src, CG_soln);
} }
Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back Rop.Dtilde(CG_soln, tmp[0]); // We actually solved Cayley preconditioned system: transform back
@ -199,7 +199,7 @@ public:
spProj(Phi, spProj_Phi, -1, Lop.Ls); spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, tmp[0], -1, 0); Lop.Omega(spProj_Phi, tmp[0], -1, 0);
G5R5(tmp[1], tmp[0]); G5R5(tmp[1], tmp[0]);
tmp[0] = zero; tmp[0] = Zero();
Solver(Lop, tmp[1], tmp[0]); Solver(Lop, tmp[1], tmp[0]);
Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back Lop.Dtilde(tmp[0], tmp[1]); // We actually solved Cayley preconditioned system: transform back
Lop.Omega(tmp[1], tmp[0], -1, 1); Lop.Omega(tmp[1], tmp[0], -1, 1);
@ -210,7 +210,7 @@ public:
spProj(Phi, spProj_Phi, 1, Rop.Ls); spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, tmp[0], 1, 0); Rop.Omega(spProj_Phi, tmp[0], 1, 0);
G5R5(tmp[1], tmp[0]); G5R5(tmp[1], tmp[0]);
tmp[0] = zero; tmp[0] = Zero();
Solver(Rop, tmp[1], tmp[0]); Solver(Rop, tmp[1], tmp[0]);
Rop.Dtilde(tmp[0], tmp[1]); Rop.Dtilde(tmp[0], tmp[1]);
Rop.Omega(tmp[1], tmp[0], 1, 1); Rop.Omega(tmp[1], tmp[0], 1, 1);
@ -238,7 +238,7 @@ public:
spProj(Phi, spProj_Phi, -1, Lop.Ls); spProj(Phi, spProj_Phi, -1, Lop.Ls);
Lop.Omega(spProj_Phi, Omega_spProj_Phi, -1, 0); Lop.Omega(spProj_Phi, Omega_spProj_Phi, -1, 0);
G5R5(CG_src, Omega_spProj_Phi); G5R5(CG_src, Omega_spProj_Phi);
spProj_Phi = zero; spProj_Phi = Zero();
Solver(Lop, CG_src, spProj_Phi); Solver(Lop, CG_src, spProj_Phi);
Lop.Dtilde(spProj_Phi, Chi); Lop.Dtilde(spProj_Phi, Chi);
G5R5(g5_R5_Chi, Chi); G5R5(g5_R5_Chi, Chi);
@ -250,7 +250,7 @@ public:
spProj(Phi, spProj_Phi, 1, Rop.Ls); spProj(Phi, spProj_Phi, 1, Rop.Ls);
Rop.Omega(spProj_Phi, Omega_spProj_Phi, 1, 0); Rop.Omega(spProj_Phi, Omega_spProj_Phi, 1, 0);
G5R5(CG_src, Omega_spProj_Phi); G5R5(CG_src, Omega_spProj_Phi);
spProj_Phi = zero; spProj_Phi = Zero();
Solver(Rop, CG_src, spProj_Phi); Solver(Rop, CG_src, spProj_Phi);
Rop.Dtilde(spProj_Phi, Chi); Rop.Dtilde(spProj_Phi, Chi);
G5R5(g5_R5_Chi, Chi); G5R5(g5_R5_Chi, Chi);

4
lib/qcd/action/pseudofermion/OneFlavourEvenOddRational.h
View File

@ -138,7 +138,7 @@ public:
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
assert(FermOp.ConstEE() == 1); assert(FermOp.ConstEE() == 1);
PhiEven = zero; PhiEven = Zero();
}; };
////////////////////////////////////////////////////// //////////////////////////////////////////////////////
@ -205,7 +205,7 @@ public:
msCG(Mpc, PhiOdd, MPhi_k); msCG(Mpc, PhiOdd, MPhi_k);
dSdU = zero; dSdU = Zero();
for (int k = 0; k < Npole; k++) { for (int k = 0; k < Npole; k++) {
RealD ak = PowerNegHalf.residues[k]; RealD ak = PowerNegHalf.residues[k];

4
lib/qcd/action/pseudofermion/OneFlavourEvenOddRationalRatio.h
View File

@ -144,7 +144,7 @@ public:
assert(NumOp.ConstEE() == 1); assert(NumOp.ConstEE() == 1);
assert(DenOp.ConstEE() == 1); assert(DenOp.ConstEE() == 1);
PhiEven = zero; PhiEven = Zero();
}; };
@ -236,7 +236,7 @@ public:
RealD ak; RealD ak;
dSdU = zero; dSdU = Zero();
// With these building blocks // With these building blocks
// //

2
lib/qcd/action/pseudofermion/OneFlavourRational.h
View File

@ -187,7 +187,7 @@ public:
msCG(MdagMOp,Phi,MPhi_k); msCG(MdagMOp,Phi,MPhi_k);
dSdU = zero; dSdU = Zero();
for(int k=0;k<Npole;k++){ for(int k=0;k<Npole;k++){
RealD ak = PowerNegHalf.residues[k]; RealD ak = PowerNegHalf.residues[k];

2
lib/qcd/action/pseudofermion/OneFlavourRationalRatio.h
View File

@ -222,7 +222,7 @@ public:
RealD ak; RealD ak;
dSdU = zero; dSdU = Zero();
// With these building blocks // With these building blocks
// //

4
lib/qcd/action/pseudofermion/TwoFlavour.h
View File

@ -111,7 +111,7 @@ public:
FermionField Y(FermOp.FermionGrid()); FermionField Y(FermOp.FermionGrid());
MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp); MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp);
X = zero; X = Zero();
ActionSolver(MdagMOp, Phi, X); ActionSolver(MdagMOp, Phi, X);
MdagMOp.Op(X, Y); MdagMOp.Op(X, Y);
@ -138,7 +138,7 @@ public:
MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp); MdagMLinearOperator<FermionOperator<Impl>, FermionField> MdagMOp(FermOp);
X = zero; X = Zero();
DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi DerivativeSolver(MdagMOp, Phi, X); // X = (MdagM)^-1 phi
MdagMOp.Op(X, Y); // Y = M X = (Mdag)^-1 phi MdagMOp.Op(X, Y); // Y = M X = (Mdag)^-1 phi

4
lib/qcd/action/pseudofermion/TwoFlavourEvenOdd.h
View File

@ -121,7 +121,7 @@ public:
SchurDifferentiableOperator<Impl> PCop(FermOp); SchurDifferentiableOperator<Impl> PCop(FermOp);
X=zero; X=Zero();
ActionSolver(PCop,PhiOdd,X); ActionSolver(PCop,PhiOdd,X);
PCop.Op(X,Y); PCop.Op(X,Y);
RealD action = norm2(Y); RealD action = norm2(Y);
@ -155,7 +155,7 @@ public:
// Our conventions really make this UdSdU; We do not differentiate wrt Udag here. // Our conventions really make this UdSdU; We do not differentiate wrt Udag here.
// So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt. // So must take dSdU - adj(dSdU) and left multiply by mom to get dS/dt.
X=zero; X=Zero();
DerivativeSolver(Mpc,PhiOdd,X); DerivativeSolver(Mpc,PhiOdd,X);
Mpc.Mpc(X,Y); Mpc.Mpc(X,Y);
Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp; Mpc.MpcDeriv(tmp , Y, X ); dSdU=tmp;

8
lib/qcd/action/pseudofermion/TwoFlavourEvenOddRatio.h
View File

@ -109,7 +109,7 @@ public:
// Odd det factors // Odd det factors
Mpc.MpcDag(etaOdd,PhiOdd); Mpc.MpcDag(etaOdd,PhiOdd);
tmp=zero; tmp=Zero();
ActionSolver(Vpc,PhiOdd,tmp); ActionSolver(Vpc,PhiOdd,tmp);
Vpc.Mpc(tmp,PhiOdd); Vpc.Mpc(tmp,PhiOdd);
@ -137,7 +137,7 @@ public:
FermionField Y(NumOp.FermionRedBlackGrid()); FermionField Y(NumOp.FermionRedBlackGrid());
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero; X=Zero();
ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi ActionSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
//Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi //Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi
// Multiply by Ydag // Multiply by Ydag
@ -145,7 +145,7 @@ public:
//RealD action = norm2(Y); //RealD action = norm2(Y);
// The EE factorised block; normally can replace with zero if det is constant (gauge field indept) // The EE factorised block; normally can replace with Zero() if det is constant (gauge field indept)
// Only really clover term that creates this. Leave the EE portion as a future to do to make most // Only really clover term that creates this. Leave the EE portion as a future to do to make most
// rapid progresss on DWF for now. // rapid progresss on DWF for now.
// //
@ -179,7 +179,7 @@ public:
//X = (Mdag M)^-1 V^dag phi //X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi //Y = (Mdag)^-1 V^dag phi
Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi Vpc.MpcDag(PhiOdd,Y); // Y= Vdag phi
X=zero; X=Zero();
DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi DerivativeSolver(Mpc,Y,X); // X= (MdagM)^-1 Vdag phi
Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi Mpc.Mpc(X,Y); // Y= Mdag^-1 Vdag phi

6
lib/qcd/action/pseudofermion/TwoFlavourRatio.h
View File

@ -94,7 +94,7 @@ public:
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp); MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(NumOp);
DenOp.Mdag(eta,Phi); // Mdag eta DenOp.Mdag(eta,Phi); // Mdag eta
tmp = zero; tmp = Zero();
ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta ActionSolver(MdagMOp,Phi,tmp); // (VdagV)^-1 Mdag eta = V^-1 Vdag^-1 Mdag eta
NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta NumOp.M(tmp,Phi); // Vdag^-1 Mdag eta
@ -116,7 +116,7 @@ public:
MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp); MdagMLinearOperator<FermionOperator<Impl> ,FermionField> MdagMOp(DenOp);
NumOp.Mdag(Phi,Y); // Y= Vdag phi NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=zero; X=Zero();
ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi ActionSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi
@ -147,7 +147,7 @@ public:
//X = (Mdag M)^-1 V^dag phi //X = (Mdag M)^-1 V^dag phi
//Y = (Mdag)^-1 V^dag phi //Y = (Mdag)^-1 V^dag phi
NumOp.Mdag(Phi,Y); // Y= Vdag phi NumOp.Mdag(Phi,Y); // Y= Vdag phi
X=zero; X=Zero();
DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi DerivativeSolver(MdagMOp,Y,X); // X= (MdagM)^-1 Vdag phi
DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi DenOp.M(X,Y); // Y= Mdag^-1 Vdag phi

2
lib/qcd/action/scalar/ScalarImpl.h
View File

@ -130,7 +130,7 @@ public:
} }
static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) { static inline void ColdConfiguration(GridParallelRNG &pRNG, Field &U) {
U = zero; U = Zero();
} }
}; };

2
lib/qcd/hmc/integrators/Integrator_algorithm.h
View File

@ -229,7 +229,7 @@ public:
Field Ufg(U.Grid()); Field Ufg(U.Grid());
Field Pfg(U.Grid()); Field Pfg(U.Grid());
Ufg = U; Ufg = U;
Pfg = zero; Pfg = Zero();
std::cout << GridLogIntegrator << "FG update " << fg_dt << " " << ep std::cout << GridLogIntegrator << "FG update " << fg_dt << " " << ep
<< std::endl; << std::endl;
// prepare_fg; no prediction/result cache for now // prepare_fg; no prediction/result cache for now

6
lib/qcd/representations/adjoint.h
View File

@ -36,7 +36,7 @@ public:
// where t^a is the generator in the fundamental // where t^a is the generator in the fundamental
// T_F is 1/2 for the fundamental representation // T_F is 1/2 for the fundamental representation
conformable(U, Uin); conformable(U, Uin);
U = zero; U = Zero();
LatticeColourMatrix tmp(Uin.Grid()); LatticeColourMatrix tmp(Uin.Grid());
Vector<typename SU<ncolour>::Matrix> ta(Dimension); Vector<typename SU<ncolour>::Matrix> ta(Dimension);
@ -76,13 +76,13 @@ public:
LatticeGaugeField RtoFundamentalProject(const LatticeField &in, LatticeGaugeField RtoFundamentalProject(const LatticeField &in,
Real scale = 1.0) const { Real scale = 1.0) const {
LatticeGaugeField out(in.Grid()); LatticeGaugeField out(in.Grid());
out = zero; out = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
LatticeColourMatrix out_mu(in.Grid()); // fundamental representation LatticeColourMatrix out_mu(in.Grid()); // fundamental representation
LatticeMatrix in_mu = peekLorentz(in, mu); LatticeMatrix in_mu = peekLorentz(in, mu);
out_mu = zero; out_mu = Zero();
typename SU<ncolour>::LatticeAlgebraVector h(in.Grid()); typename SU<ncolour>::LatticeAlgebraVector h(in.Grid());
projectOnAlgebra(h, in_mu, double(Nc) * 2.0); // factor C(r)/C(fund) projectOnAlgebra(h, in_mu, double(Nc) * 2.0); // factor C(r)/C(fund)

6
lib/qcd/representations/two_index.h
View File

@ -39,7 +39,7 @@ public:
// get the U in TwoIndexRep // get the U in TwoIndexRep
// (U)_{(ij)(lk)} = tr [ adj(e^(ij)) U e^(lk) transpose(U) ] // (U)_{(ij)(lk)} = tr [ adj(e^(ij)) U e^(lk) transpose(U) ]
conformable(U, Uin); conformable(U, Uin);
U = zero; U = Zero();
LatticeColourMatrix tmp(Uin.Grid()); LatticeColourMatrix tmp(Uin.Grid());
Vector<typename SU<ncolour>::Matrix> eij(Dimension); Vector<typename SU<ncolour>::Matrix> eij(Dimension);
@ -62,13 +62,13 @@ public:
LatticeGaugeField RtoFundamentalProject(const LatticeField &in, LatticeGaugeField RtoFundamentalProject(const LatticeField &in,
Real scale = 1.0) const { Real scale = 1.0) const {
LatticeGaugeField out(in.Grid()); LatticeGaugeField out(in.Grid());
out = zero; out = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
LatticeColourMatrix out_mu(in.Grid()); // fundamental representation LatticeColourMatrix out_mu(in.Grid()); // fundamental representation
LatticeMatrix in_mu = peekLorentz(in, mu); LatticeMatrix in_mu = peekLorentz(in, mu);
out_mu = zero; out_mu = Zero();
typename SU<ncolour>::LatticeAlgebraVector h(in.Grid()); typename SU<ncolour>::LatticeAlgebraVector h(in.Grid());
projectOnAlgebra(h, in_mu, double(Nc + 2 * S)); // factor T(r)/T(fund) projectOnAlgebra(h, in_mu, double(Nc + 2 * S)); // factor T(r)/T(fund)

6
lib/qcd/smearing/APEsmearing.h
View File

@ -66,10 +66,10 @@ public:
GridBase *grid = U.Grid(); GridBase *grid = U.Grid();
GaugeLinkField Cup(grid), tmp_stpl(grid); GaugeLinkField Cup(grid), tmp_stpl(grid);
WilsonLoops<Gimpl> WL; WilsonLoops<Gimpl> WL;
u_smr = zero; u_smr = Zero();
for(int mu=0; mu<Nd; ++mu){ for(int mu=0; mu<Nd; ++mu){
Cup = zero; Cup = Zero();
for(int nu=0; nu<Nd; ++nu){ for(int nu=0; nu<Nd; ++nu){
if (nu != mu) { if (nu != mu) {
// get the staple in direction mu, nu // get the staple in direction mu, nu
@ -132,7 +132,7 @@ public:
//-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x) //-r_munu*U_nu(x+mu)*Udag_mu(x+nu)*Lambda_mu(x+nu)*Udag_nu(x)
Gimpl::AddLink(SigmaTerm, temp_Sigma, mu); Gimpl::AddLink(SigmaTerm, temp_Sigma, mu);
staple = zero; staple = Zero();
sh_field = Cshift(U_nu, mu, 1); sh_field = Cshift(U_nu, mu, 1);
temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu; temp_Sigma = -rho_munu*adj(sh_field)*adj(U_mu)*iLambda_mu*U_nu;

6
lib/qcd/smearing/GaugeConfiguration.h
View File

@ -64,7 +64,7 @@ private:
if (smearingLevels > 0) { if (smearingLevels > 0) {
std::cout << GridLogDebug std::cout << GridLogDebug
<< "[SmearedConfiguration] Filling SmearedSet\n"; << "[SmearedConfiguration] Filling SmearedSet\n";
GaugeField previous_u(ThinLinks->_grid); GaugeField previous_u(ThinLinks->Grid());
previous_u = *ThinLinks; previous_u = *ThinLinks;
for (int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl) { for (int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl) {
@ -89,8 +89,8 @@ private:
GaugeLinkField GaugeKmu(grid), Cmu(grid); GaugeLinkField GaugeKmu(grid), Cmu(grid);
StoutSmearing.BaseSmear(C, GaugeK); StoutSmearing.BaseSmear(C, GaugeK);
SigmaK = zero; SigmaK = Zero();
iLambda = zero; iLambda = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
Cmu = peekLorentz(C, mu); Cmu = peekLorentz(C, mu);

16
lib/qcd/spin/TwoSpinor.h
View File

@ -159,15 +159,15 @@ template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spProj
//typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE; //typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE;
rfspin(0)=fspin(0); rfspin(0)=fspin(0);
rfspin(1)=fspin(1); rfspin(1)=fspin(1);
rfspin(2)=zero; rfspin(2)=Zero();
rfspin(3)=zero; rfspin(3)=Zero();
} }
// template<class vtype> strong_inline void fspProj5m (iVector<vtype,Ns> &rfspin,const iVector<vtype,Ns> &fspin) // template<class vtype> strong_inline void fspProj5m (iVector<vtype,Ns> &rfspin,const iVector<vtype,Ns> &fspin)
template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spProj5m (iVector<vtype,Ns> &rfspin,const iVector<vtype,Ns> &fspin) template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spProj5m (iVector<vtype,Ns> &rfspin,const iVector<vtype,Ns> &fspin)
{ {
//typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE; //typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE;
rfspin(0)=zero; rfspin(0)=Zero();
rfspin(1)=zero; rfspin(1)=Zero();
rfspin(2)=fspin(2); rfspin(2)=fspin(2);
rfspin(3)=fspin(3); rfspin(3)=fspin(3);
} }
@ -339,14 +339,14 @@ template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spReco
//typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE; //typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE;
fspin(0)=hspin(0)+hspin(0); // add is lower latency than mul fspin(0)=hspin(0)+hspin(0); // add is lower latency than mul
fspin(1)=hspin(1)+hspin(1); // probably no measurable diffence though fspin(1)=hspin(1)+hspin(1); // probably no measurable diffence though
fspin(2)=zero; fspin(2)=Zero();
fspin(3)=zero; fspin(3)=Zero();
} }
template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spRecon5m (iVector<vtype,Ns> &fspin,const iVector<vtype,Nhs> &hspin) template<class vtype,IfSpinor<iVector<vtype,Ns> > = 0> strong_inline void spRecon5m (iVector<vtype,Ns> &fspin,const iVector<vtype,Nhs> &hspin)
{ {
//typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE; //typename std::enable_if<matchGridTensorIndex<iVector<vtype,Ns>,SpinorIndex>::value,iVector<vtype,Ns> >::type *SFINAE;
fspin(0)=zero; fspin(0)=Zero();
fspin(1)=zero; fspin(1)=Zero();
fspin(2)=hspin(0)+hspin(0); fspin(2)=hspin(0)+hspin(0);
fspin(3)=hspin(1)+hspin(1); fspin(3)=hspin(1)+hspin(1);
} }

6
lib/qcd/utils/CovariantLaplacian.h
View File

@ -113,7 +113,7 @@ public:
GaugeLinkField sum(in.Grid()); GaugeLinkField sum(in.Grid());
for (int nu = 0; nu < Nd; nu++) { for (int nu = 0; nu < Nd; nu++) {
sum = zero; sum = Zero();
GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu); GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
GaugeLinkField out_nu(out.Grid()); GaugeLinkField out_nu(out.Grid());
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
@ -132,7 +132,7 @@ public:
for (int mu = 0; mu < Nd; mu++){ for (int mu = 0; mu < Nd; mu++){
GaugeLinkField der_mu(der.Grid()); GaugeLinkField der_mu(der.Grid());
der_mu = zero; der_mu = Zero();
for (int nu = 0; nu < Nd; nu++){ for (int nu = 0; nu < Nd; nu++){
GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu); GaugeLinkField in_nu = PeekIndex<LorentzIndex>(in, nu);
der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu; der_mu += U[mu] * Cshift(in_nu, mu, 1) * adj(U[mu]) * in_nu;
@ -151,7 +151,7 @@ public:
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
GaugeLinkField der_mu(der.Grid()); GaugeLinkField der_mu(der.Grid());
der_mu = zero; der_mu = Zero();
for (int nu = 0; nu < Nd; nu++) { for (int nu = 0; nu < Nd; nu++) {
GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu); GaugeLinkField left_nu = PeekIndex<LorentzIndex>(left, nu);
GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu); GaugeLinkField right_nu = PeekIndex<LorentzIndex>(right, nu);

4
lib/qcd/utils/GaugeFix.h
View File

@ -46,7 +46,7 @@ public:
} }
} }
static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu) { static void DmuAmu(const std::vector<GaugeMat> &A,GaugeMat &dmuAmu) {
dmuAmu=zero; dmuAmu=Zero();
for(int mu=0;mu<Nd;mu++){ for(int mu=0;mu<Nd;mu++){
dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1); dmuAmu = dmuAmu + A[mu] - Cshift(A[mu],mu,-1);
} }
@ -123,7 +123,7 @@ public:
FFT theFFT((GridCartesian *)grid); FFT theFFT((GridCartesian *)grid);
LatticeComplex Fp(grid); LatticeComplex Fp(grid);
LatticeComplex psq(grid); psq=zero; LatticeComplex psq(grid); psq=Zero();
LatticeComplex pmu(grid); LatticeComplex pmu(grid);
LatticeComplex one(grid); one = Complex(1.0,0.0); LatticeComplex one(grid); one = Complex(1.0,0.0);

16
lib/qcd/utils/Metric.h
View File

@ -63,10 +63,10 @@ public:
// do nothing // do nothing
} }
virtual void MDeriv(const Field& in, Field& out){ virtual void MDeriv(const Field& in, Field& out){
out = zero; out = Zero();
} }
virtual void MDeriv(const Field& left, const Field& right, Field& out){ virtual void MDeriv(const Field& left, const Field& right, Field& out){
out = zero; out = Zero();
} }
}; };
@ -119,10 +119,10 @@ public:
// Correct // Correct
RealD MomentaAction(){ RealD MomentaAction(){
MomentaField inv(Mom.Grid()); MomentaField inv(Mom.Grid());
inv = zero; inv = Zero();
M.Minv(Mom, inv); M.Minv(Mom, inv);
LatticeComplex Hloc(Mom.Grid()); LatticeComplex Hloc(Mom.Grid());
Hloc = zero; Hloc = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
// This is not very general // This is not very general
// hide in the metric // hide in the metric
@ -157,7 +157,7 @@ public:
// with respect to the gauge field // with respect to the gauge field
MomentaField MDer(in.Grid()); MomentaField MDer(in.Grid());
MomentaField X(in.Grid()); MomentaField X(in.Grid());
X = zero; X = Zero();
M.Minv(in, X); // X = G in M.Minv(in, X); // X = G in
M.MDeriv(X, MDer); // MDer = U * dS/dU M.MDeriv(X, MDer); // MDer = U * dS/dU
der = Implementation::projectForce(MDer); // Ta if gauge fields der = Implementation::projectForce(MDer); // Ta if gauge fields
@ -165,12 +165,12 @@ public:
} }
void AuxiliaryFieldsDerivative(MomentaField& der){ void AuxiliaryFieldsDerivative(MomentaField& der){
der = zero; der = Zero();
if (1){ if (1){
// Auxiliary fields // Auxiliary fields
MomentaField der_temp(der.Grid()); MomentaField der_temp(der.Grid());
MomentaField X(der.Grid()); MomentaField X(der.Grid());
X=zero; X=Zero();
//M.M(AuxMom, X); // X = M Aux //M.M(AuxMom, X); // X = M Aux
// Two derivative terms // Two derivative terms
// the Mderiv need separation of left and right terms // the Mderiv need separation of left and right terms
@ -185,7 +185,7 @@ public:
} }
void DerivativeP(MomentaField& der){ void DerivativeP(MomentaField& der){
der = zero; der = Zero();
M.Minv(Mom, der); M.Minv(Mom, der);
// is the projection necessary here? // is the projection necessary here?
// no for fields in the algebra // no for fields in the algebra

30
lib/qcd/utils/SUn.h
View File

@ -162,7 +162,7 @@ public:
template <class cplx> template <class cplx>
static void generatorSigmaY(int su2Index, iSUnMatrix<cplx> &ta) { static void generatorSigmaY(int su2Index, iSUnMatrix<cplx> &ta) {
ta = zero; ta = Zero();
int i1, i2; int i1, i2;
su2SubGroupIndex(i1, i2, su2Index); su2SubGroupIndex(i1, i2, su2Index);
ta()()(i1, i2) = 1.0; ta()()(i1, i2) = 1.0;
@ -172,7 +172,7 @@ public:
template <class cplx> template <class cplx>
static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) { static void generatorSigmaX(int su2Index, iSUnMatrix<cplx> &ta) {
ta = zero; ta = Zero();
cplx i(0.0, 1.0); cplx i(0.0, 1.0);
int i1, i2; int i1, i2;
su2SubGroupIndex(i1, i2, su2Index); su2SubGroupIndex(i1, i2, su2Index);
@ -184,7 +184,7 @@ public:
template <class cplx> template <class cplx>
static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) { static void generatorDiagonal(int diagIndex, iSUnMatrix<cplx> &ta) {
// diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...)
ta = zero; ta = Zero();
int k = diagIndex + 1; // diagIndex starts from 0 int k = diagIndex + 1; // diagIndex starts from 0
for (int i = 0; i <= diagIndex; i++) { // k iterations for (int i = 0; i <= diagIndex; i++) { // k iterations
ta()()(i, i) = 1.0; ta()()(i, i) = 1.0;
@ -299,11 +299,11 @@ public:
LatticeComplex ones(grid); LatticeComplex ones(grid);
ones = 1.0; ones = 1.0;
LatticeComplex zeros(grid); LatticeComplex zeros(grid);
zeros = zero; zeros = Zero();
LatticeReal rones(grid); LatticeReal rones(grid);
rones = 1.0; rones = 1.0;
LatticeReal rzeros(grid); LatticeReal rzeros(grid);
rzeros = zero; rzeros = Zero();
LatticeComplex udet(grid); // determinant of real(staple) LatticeComplex udet(grid); // determinant of real(staple)
LatticeInteger mask_true(grid); LatticeInteger mask_true(grid);
mask_true = 1; mask_true = 1;
@ -432,13 +432,13 @@ public:
RealD numSites = sum(rtmp); RealD numSites = sum(rtmp);
RealD numAccepted; RealD numAccepted;
LatticeInteger Accepted(grid); LatticeInteger Accepted(grid);
Accepted = zero; Accepted = Zero();
LatticeInteger newlyAccepted(grid); LatticeInteger newlyAccepted(grid);
std::vector<LatticeReal> xr(4, grid); std::vector<LatticeReal> xr(4, grid);
std::vector<LatticeReal> a(4, grid); std::vector<LatticeReal> a(4, grid);
LatticeReal d(grid); LatticeReal d(grid);
d = zero; d = Zero();
LatticeReal alpha(grid); LatticeReal alpha(grid);
// std::cout<<GridLogMessage<<"xi "<<xi <<std::endl; // std::cout<<GridLogMessage<<"xi "<<xi <<std::endl;
@ -475,7 +475,7 @@ public:
LatticeInteger ione(grid); LatticeInteger ione(grid);
ione = 1; ione = 1;
LatticeInteger izero(grid); LatticeInteger izero(grid);
izero = zero; izero = Zero();
newlyAccepted = where(xrsq < thresh, ione, izero); newlyAccepted = where(xrsq < thresh, ione, izero);
Accepted = where(newlyAccepted, newlyAccepted, Accepted); Accepted = where(newlyAccepted, newlyAccepted, Accepted);
@ -490,7 +490,7 @@ public:
} while ((numAccepted < numSites) && (hit < nheatbath)); } while ((numAccepted < numSites) && (hit < nheatbath));
// G. Set a0 = 1 - d; // G. Set a0 = 1 - d;
a[0] = zero; a[0] = Zero();
a[0] = where(wheremask, 1.0 - d, a[0]); a[0] = where(wheremask, 1.0 - d, a[0]);
////////////////////////////////////////// //////////////////////////////////////////
@ -528,7 +528,7 @@ public:
// Debug Checks // Debug Checks
// SU2 check // SU2 check
LatticeSU2Matrix check(grid); // rotated matrix after hb LatticeSU2Matrix check(grid); // rotated matrix after hb
u = zero; u = Zero();
check = ua * adj(ua) - 1.0; check = ua * adj(ua) - 1.0;
check = where(Accepted, check, u); check = where(Accepted, check, u);
assert(norm2(check) < 1.0e-4); assert(norm2(check) < 1.0e-4);
@ -538,7 +538,7 @@ public:
assert(norm2(check) < 1.0e-4); assert(norm2(check) < 1.0e-4);
LatticeMatrix Vcheck(grid); LatticeMatrix Vcheck(grid);
Vcheck = zero; Vcheck = Zero();
Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck);
// std::cout<<GridLogMessage << "SU3 check " <<norm2(Vcheck)<<std::endl; // std::cout<<GridLogMessage << "SU3 check " <<norm2(Vcheck)<<std::endl;
assert(norm2(Vcheck) < 1.0e-4); assert(norm2(Vcheck) < 1.0e-4);
@ -622,7 +622,7 @@ public:
ComplexD cone(1.0, 0.0); ComplexD cone(1.0, 0.0);
MatrixType ta; MatrixType ta;
lie = zero; lie = Zero();
for (int a = 0; a < AdjointDimension; a++) { for (int a = 0; a < AdjointDimension; a++) {
random(pRNG, ca); random(pRNG, ca);
@ -647,7 +647,7 @@ public:
Complex ci(0.0, scale); Complex ci(0.0, scale);
Matrix ta; Matrix ta;
out = zero; out = Zero();
for (int a = 0; a < AdjointDimension; a++) { for (int a = 0; a < AdjointDimension; a++) {
gaussian(pRNG, ca); gaussian(pRNG, ca);
generator(a, ta); generator(a, ta);
@ -665,7 +665,7 @@ public:
LatticeMatrix la(grid); LatticeMatrix la(grid);
Matrix ta; Matrix ta;
out = zero; out = Zero();
for (int a = 0; a < AdjointDimension; a++) { for (int a = 0; a < AdjointDimension; a++) {
generator(a, ta); generator(a, ta);
la = peekColour(h, a) * timesI(ta) * scale; la = peekColour(h, a) * timesI(ta) * scale;
@ -708,7 +708,7 @@ public:
// inverse operation: FundamentalLieAlgebraMatrix // inverse operation: FundamentalLieAlgebraMatrix
static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) { static void projectOnAlgebra(LatticeAlgebraVector &h_out, const LatticeMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
h_out = zero; h_out = Zero();
Matrix Ta; Matrix Ta;
for (int a = 0; a < AdjointDimension; a++) { for (int a = 0; a < AdjointDimension; a++) {

8
lib/qcd/utils/SUnAdjoint.h
View File

@ -61,7 +61,7 @@ public:
static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) { static void generator(int Index, iSUnAdjointMatrix<cplx> &iAdjTa) {
// returns i(T_Adj)^index necessary for the projectors // returns i(T_Adj)^index necessary for the projectors
// see definitions above // see definitions above
iAdjTa = zero; iAdjTa = Zero();
Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta(ncolour * ncolour - 1); Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > ta(ncolour * ncolour - 1);
typename SU<ncolour>::template iSUnMatrix<cplx> tmp; typename SU<ncolour>::template iSUnMatrix<cplx> tmp;
@ -118,7 +118,7 @@ public:
LatticeAdjMatrix la(grid); LatticeAdjMatrix la(grid);
AMatrix iTa; AMatrix iTa;
out = zero; out = Zero();
for (int a = 0; a < Dimension; a++) { for (int a = 0; a < Dimension; a++) {
generator(a, iTa); generator(a, iTa);
la = peekColour(h, a) * iTa; la = peekColour(h, a) * iTa;
@ -130,7 +130,7 @@ public:
// Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 )
static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) { static void projectOnAlgebra(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
h_out = zero; h_out = Zero();
AMatrix iTa; AMatrix iTa;
Real coefficient = - 1.0/(ncolour) * scale;// 1/Nc for the normalization of the trace in the adj rep Real coefficient = - 1.0/(ncolour) * scale;// 1/Nc for the normalization of the trace in the adj rep
@ -145,7 +145,7 @@ public:
static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) { static void projector(typename SU<ncolour>::LatticeAlgebraVector &h_out, const LatticeAdjMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
static std::vector<AMatrix> iTa(Dimension); // to store the generators static std::vector<AMatrix> iTa(Dimension); // to store the generators
h_out = zero; h_out = Zero();
static bool precalculated = false; static bool precalculated = false;
if (!precalculated){ if (!precalculated){
precalculated = true; precalculated = true;

14
lib/qcd/utils/SUnTwoIndex.h
View File

@ -71,7 +71,7 @@ public:
static void base(int Index, iSUnMatrix<cplx> &eij) { static void base(int Index, iSUnMatrix<cplx> &eij) {
// returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R // returns (e)^(ij)_{kl} necessary for change of base U_F -> U_R
assert(Index < NumGenerators); assert(Index < NumGenerators);
eij = zero; eij = Zero();
// for the linearisation of the 2 indexes // for the linearisation of the 2 indexes
static int a[ncolour * (ncolour - 1) / 2][2]; // store the a <-> i,j static int a[ncolour * (ncolour - 1) / 2][2]; // store the a <-> i,j
@ -97,14 +97,14 @@ public:
template <class cplx> template <class cplx>
static void baseDiagonal(int Index, iSUnMatrix<cplx> &eij) { static void baseDiagonal(int Index, iSUnMatrix<cplx> &eij) {
eij = zero; eij = Zero();
eij()()(Index - ncolour * (ncolour - 1) / 2, eij()()(Index - ncolour * (ncolour - 1) / 2,
Index - ncolour * (ncolour - 1) / 2) = 1.0; Index - ncolour * (ncolour - 1) / 2) = 1.0;
} }
template <class cplx> template <class cplx>
static void baseOffDiagonal(int i, int j, iSUnMatrix<cplx> &eij) { static void baseOffDiagonal(int i, int j, iSUnMatrix<cplx> &eij) {
eij = zero; eij = Zero();
for (int k = 0; k < ncolour; k++) for (int k = 0; k < ncolour; k++)
for (int l = 0; l < ncolour; l++) for (int l = 0; l < ncolour; l++)
eij()()(l, k) = delta(i, k) * delta(j, l) + eij()()(l, k) = delta(i, k) * delta(j, l) +
@ -130,7 +130,7 @@ public:
ncolour * ncolour - 1); ncolour * ncolour - 1);
Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > eij(Dimension); Vector<typename SU<ncolour>::template iSUnMatrix<cplx> > eij(Dimension);
typename SU<ncolour>::template iSUnMatrix<cplx> tmp; typename SU<ncolour>::template iSUnMatrix<cplx> tmp;
i2indTa = zero; i2indTa = Zero();
for (int a = 0; a < ncolour * ncolour - 1; a++) for (int a = 0; a < ncolour * ncolour - 1; a++)
SU<ncolour>::generator(a, ta[a]); SU<ncolour>::generator(a, ta[a]);
@ -203,7 +203,7 @@ public:
LatticeTwoIndexMatrix la(grid); LatticeTwoIndexMatrix la(grid);
TIMatrix i2indTa; TIMatrix i2indTa;
out = zero; out = Zero();
for (int a = 0; a < ncolour * ncolour - 1; a++) { for (int a = 0; a < ncolour * ncolour - 1; a++) {
generator(a, i2indTa); generator(a, i2indTa);
la = peekColour(h, a) * i2indTa; la = peekColour(h, a) * i2indTa;
@ -218,7 +218,7 @@ public:
typename SU<ncolour>::LatticeAlgebraVector &h_out, typename SU<ncolour>::LatticeAlgebraVector &h_out,
const LatticeTwoIndexMatrix &in, Real scale = 1.0) { const LatticeTwoIndexMatrix &in, Real scale = 1.0) {
conformable(h_out, in); conformable(h_out, in);
h_out = zero; h_out = Zero();
TIMatrix i2indTa; TIMatrix i2indTa;
Real coefficient = -2.0 / (ncolour + 2 * S) * scale; Real coefficient = -2.0 / (ncolour + 2 * S) * scale;
// 2/(Nc +/- 2) for the normalization of the trace in the two index rep // 2/(Nc +/- 2) for the normalization of the trace in the two index rep
@ -236,7 +236,7 @@ public:
conformable(h_out, in); conformable(h_out, in);
// to store the generators // to store the generators
static std::vector<TIMatrix> i2indTa(ncolour * ncolour -1); static std::vector<TIMatrix> i2indTa(ncolour * ncolour -1);
h_out = zero; h_out = Zero();
static bool precalculated = false; static bool precalculated = false;
if (!precalculated) { if (!precalculated) {
precalculated = true; precalculated = true;

16
lib/qcd/utils/WilsonLoops.h
View File

@ -86,7 +86,7 @@ public:
static void sitePlaquette(ComplexField &Plaq, static void sitePlaquette(ComplexField &Plaq,
const std::vector<GaugeMat> &U) { const std::vector<GaugeMat> &U) {
ComplexField sitePlaq(U[0].Grid()); ComplexField sitePlaq(U[0].Grid());
Plaq = zero; Plaq = Zero();
for (int mu = 1; mu < Nd; mu++) { for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) { for (int nu = 0; nu < mu; nu++) {
traceDirPlaquette(sitePlaq, U, mu, nu); traceDirPlaquette(sitePlaq, U, mu, nu);
@ -130,7 +130,7 @@ public:
std::vector<GaugeMat> U(Nd, Umu.Grid()); std::vector<GaugeMat> U(Nd, Umu.Grid());
ComplexField Tr(Umu.Grid()); ComplexField Tr(Umu.Grid());
Tr = zero; Tr = Zero();
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu); U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
Tr = Tr + trace(U[mu]); Tr = Tr + trace(U[mu]);
@ -156,7 +156,7 @@ public:
for (int d = 0; d < Nd; d++) { for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d); U[d] = PeekIndex<LorentzIndex>(Umu, d);
} }
staple = zero; staple = Zero();
if (nu != mu) { if (nu != mu) {
@ -197,7 +197,7 @@ public:
// this operation is taking too much time // this operation is taking too much time
U[d] = PeekIndex<LorentzIndex>(Umu, d); U[d] = PeekIndex<LorentzIndex>(Umu, d);
} }
staple = zero; staple = Zero();
GaugeMat tmp1(grid); GaugeMat tmp1(grid);
GaugeMat tmp2(grid); GaugeMat tmp2(grid);
@ -225,7 +225,7 @@ public:
for (int d = 0; d < Nd; d++) { for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d); U[d] = PeekIndex<LorentzIndex>(Umu, d);
} }
staple = zero; staple = Zero();
for (int nu = 0; nu < Nd; nu++) { for (int nu = 0; nu < Nd; nu++) {
@ -385,7 +385,7 @@ public:
static void siteRectangle(ComplexField &Rect, static void siteRectangle(ComplexField &Rect,
const std::vector<GaugeMat> &U) { const std::vector<GaugeMat> &U) {
ComplexField siteRect(U[0].Grid()); ComplexField siteRect(U[0].Grid());
Rect = zero; Rect = Zero();
for (int mu = 1; mu < Nd; mu++) { for (int mu = 1; mu < Nd; mu++) {
for (int nu = 0; nu < mu; nu++) { for (int nu = 0; nu < mu; nu++) {
traceDirRectangle(siteRect, U, mu, nu); traceDirRectangle(siteRect, U, mu, nu);
@ -443,7 +443,7 @@ public:
static void RectStapleOptimised(GaugeMat &Stap, std::vector<GaugeMat> &U2, static void RectStapleOptimised(GaugeMat &Stap, std::vector<GaugeMat> &U2,
std::vector<GaugeMat> &U, int mu) { std::vector<GaugeMat> &U, int mu) {
Stap = zero; Stap = Zero();
GridBase *grid = U[0].Grid(); GridBase *grid = U[0].Grid();
@ -529,7 +529,7 @@ public:
U[d] = PeekIndex<LorentzIndex>(Umu, d); U[d] = PeekIndex<LorentzIndex>(Umu, d);
} }
Stap = zero; Stap = Zero();
for (int nu = 0; nu < Nd; nu++) { for (int nu = 0; nu < Nd; nu++) {
if (nu != mu) { if (nu != mu) {