mirror of
https://github.com/paboyle/Grid.git
synced 2024-11-10 07:55:35 +00:00
Merge branch 'develop' of https://github.com/paboyle/Grid into develop
This commit is contained in:
commit
6957b0b58a
@ -354,6 +354,6 @@ public:
|
||||
}
|
||||
};
|
||||
|
||||
NAMESPACE_END(QCD);
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
#endif
|
||||
|
@ -10,6 +10,7 @@ Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Gianluca Filaci <g.filaci@ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
@ -54,6 +55,10 @@ CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
auto chi = chi_i.View();
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
|
||||
int Ls =this->Ls;
|
||||
|
||||
// 10 = 3 complex mult + 2 complex add
|
||||
@ -71,7 +76,7 @@ CayleyFermion5D<Impl>::M5D(const FermionField &psi_i,
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5m(tmp1,psi(idx_u));
|
||||
spProj5p(tmp2,psi(idx_l));
|
||||
coalescedWrite(chi[ss+s],diag[s]*phi(ss+s)+upper[s]*tmp1+lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s],pdiag[s]*phi(ss+s)+pupper[s]*tmp1+plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
M5Dtime+=usecond();
|
||||
@ -93,6 +98,10 @@ CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
|
||||
auto chi = chi_i.View();
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
|
||||
int Ls=this->Ls;
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
@ -109,7 +118,7 @@ CayleyFermion5D<Impl>::M5Ddag(const FermionField &psi_i,
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5p(tmp1,psi(idx_u));
|
||||
spProj5m(tmp2,psi(idx_l));
|
||||
coalescedWrite(chi[ss+s],diag[s]*phi(ss+s)+upper[s]*tmp1+lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s],pdiag[s]*phi(ss+s)+pupper[s]*tmp1+plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
M5Dtime+=usecond();
|
||||
@ -139,34 +148,36 @@ CayleyFermion5D<Impl>::MooeeInv (const FermionField &psi_i, FermionField &chi
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// flops = 12*2*Ls + 12*2*Ls + 3*12*Ls + 12*2*Ls = 12*Ls * (9) = 108*Ls flops
|
||||
// Apply (L^{\prime})^{-1}
|
||||
coalescedWrite(chi[ss],psi(ss)); // chi[0]=psi[0]
|
||||
for(int s=1;s<Ls;s++){
|
||||
spProj5p(tmp,chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s] , psi(ss+s)-plee[s-1]*tmp);
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (L^{\prime})^{-1} L_m^{-1}
|
||||
res = psi(ss);
|
||||
spProj5m(tmp,res);
|
||||
acc = pleem[0]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= plee[s-1]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
acc += pleem[s]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1) - plee[Ls-2]*tmp - acc;
|
||||
|
||||
// L_m^{-1}
|
||||
for (int s=0;s<Ls-1;s++){ // Chi[ee] = 1 - sum[s<Ls-1] -pleem[s]P_- chi
|
||||
spProj5m(tmp,chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pleem[s]*tmp);
|
||||
}
|
||||
|
||||
// U_m^{-1} D^{-1}
|
||||
for (int s=0;s<Ls-1;s++){
|
||||
// Chi[s] + 1/d chi[s]
|
||||
spProj5p(tmp,chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s)-(pueem[s]/pdee[Ls-1])*tmp);
|
||||
}
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
|
||||
// Apply U^{-1}
|
||||
// Apply U_m^{-1} D^{-1} U^{-1}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
coalescedWrite(chi[ss+Ls-1],res);
|
||||
spProj5p(acc,res);
|
||||
spProj5m(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
spProj5m(tmp,chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - puee[s]*tmp);
|
||||
res = (1.0/pdee[s])*chi(ss+s) - puee[s]*tmp - pueem[s]*acc;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
@ -201,31 +212,36 @@ CayleyFermion5D<Impl>::MooeeInvDag (const FermionField &psi_i, FermionField &chi
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// Apply (U^{\prime})^{-dagger}
|
||||
coalescedWrite(chi[ss],psi(ss));
|
||||
for (int s=1;s<Ls;s++){
|
||||
spProj5m(tmp,chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s)-conjugate(puee[s-1])*tmp);
|
||||
}
|
||||
// U_m^{-\dagger}
|
||||
for (int s=0;s<Ls-1;s++){
|
||||
spProj5p(tmp,chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - conjugate(pueem[s])*tmp);
|
||||
}
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (U^{\prime})^{-dagger} U_m^{-\dagger}
|
||||
res = psi(ss);
|
||||
spProj5p(tmp,res);
|
||||
acc = conjugate(pueem[0])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
// L_m^{-\dagger} D^{-dagger}
|
||||
for (int s=0;s<Ls-1;s++){
|
||||
spProj5m(tmp,chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], conjugate(1.0/pdee[s])*chi(ss+s)-conjugate(pleem[s]/pdee[Ls-1])*tmp);
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= conjugate(puee[s-1])*tmp;
|
||||
spProj5p(tmp,res);
|
||||
acc += conjugate(pueem[s])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
coalescedWrite(chi[ss+Ls-1], conjugate(1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
res = psi(ss+Ls-1) - conjugate(puee[Ls-2])*tmp - acc;
|
||||
|
||||
// Apply L^{-dagger}
|
||||
// Apply L_m^{-\dagger} D^{-dagger} L^{-dagger}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
coalescedWrite(chi[ss+Ls-1],res);
|
||||
spProj5m(acc,res);
|
||||
spProj5p(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
spProj5p(tmp,chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - conjugate(plee[s])*tmp);
|
||||
res = (1.0/pdee[s])*chi(ss+s) - conjugate(plee[s])*tmp - conjugate(pleem[s])*acc;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
MooeeInvTime+=usecond();
|
||||
|
@ -11,6 +11,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: David Murphy <dmurphy@phys.columbia.edu>
|
||||
Author: Gianluca Filaci <g.filaci@ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
@ -49,6 +50,9 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
|
||||
auto psi = psi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
this->M5Dtime -= usecond();
|
||||
@ -63,7 +67,7 @@ void DomainWallEOFAFermion<Impl>::M5D(const FermionField& psi_i, const FermionFi
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5m(tmp1, psi(idx_u));
|
||||
spProj5p(tmp2, psi(idx_l));
|
||||
coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
|
||||
@ -82,6 +86,9 @@ void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const Fermio
|
||||
auto phi = phi_i.View();
|
||||
auto chi = chi_i.View();
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
@ -97,7 +104,7 @@ void DomainWallEOFAFermion<Impl>::M5Ddag(const FermionField& psi_i, const Fermio
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5p(tmp1, psi(idx_u));
|
||||
spProj5m(tmp2, psi(idx_l));
|
||||
coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
|
||||
@ -124,36 +131,37 @@ void DomainWallEOFAFermion<Impl>::MooeeInv(const FermionField& psi_i, FermionFie
|
||||
this->MooeeInvTime -= usecond();
|
||||
uint64_t nloop=grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
auto ss=sss*Ls;
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp1,tmp2;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// flops = 12*2*Ls + 12*2*Ls + 3*12*Ls + 12*2*Ls = 12*Ls * (9) = 108*Ls flops
|
||||
// Apply (L^{\prime})^{-1}
|
||||
coalescedWrite(chi[ss],psi(ss)); // chi[0]=psi[0]
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5p(tmp1, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s) - plee[s-1]*tmp1);
|
||||
// Apply (L^{\prime})^{-1} L_m^{-1}
|
||||
res = psi(ss);
|
||||
spProj5m(tmp,res);
|
||||
acc = pleem[0]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= plee[s-1]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
acc += pleem[s]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1) - plee[Ls-2]*tmp - acc;
|
||||
|
||||
// L_m^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
|
||||
spProj5m(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pleem[s]*tmp1);
|
||||
}
|
||||
|
||||
// U_m^{-1} D^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[s] + 1/d chi[s]
|
||||
spProj5p(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s) - (pueem[s]/pdee[Ls])*tmp1);
|
||||
}
|
||||
spProj5m(tmp2, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+Ls-1],(1.0/pdee[Ls])*tmp1 + (1.0/pdee[Ls-1])*tmp2);
|
||||
|
||||
// Apply U^{-1}
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
spProj5m(tmp1, chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - puee[s]*tmp1);
|
||||
// Apply U_m^{-1} D^{-1} U^{-1}
|
||||
acc = (1.0/pdee[Ls ])*res;
|
||||
tmp = (1.0/pdee[Ls-1])*res;
|
||||
spProj5p(acc,acc);
|
||||
spProj5m(tmp,tmp);
|
||||
coalescedWrite(chi[ss+Ls-1], acc + tmp);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = (1.0/pdee[s])*chi(ss+s) - puee[s]*tmp - pueem[s]*acc;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
this->MooeeInvTime += usecond();
|
||||
@ -168,56 +176,50 @@ void DomainWallEOFAFermion<Impl>::MooeeInvDag(const FermionField& psi_i, Fermion
|
||||
auto chi = chi_i.View();
|
||||
int Ls = this->Ls;
|
||||
|
||||
auto plee = & this->lee[0];
|
||||
auto pdee = & this->dee[0];
|
||||
auto puee = & this->uee[0];
|
||||
|
||||
auto pleem = & this->leem[0];
|
||||
auto pueem = & this->ueem[0];
|
||||
|
||||
assert(psi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
Vector<Coeff_t> ueec(Ls);
|
||||
Vector<Coeff_t> deec(Ls+1);
|
||||
Vector<Coeff_t> leec(Ls);
|
||||
Vector<Coeff_t> ueemc(Ls);
|
||||
Vector<Coeff_t> leemc(Ls);
|
||||
|
||||
for(int s=0; s<ueec.size(); s++){
|
||||
ueec[s] = conjugate(this->uee[s]);
|
||||
deec[s] = conjugate(this->dee[s]);
|
||||
leec[s] = conjugate(this->lee[s]);
|
||||
ueemc[s] = conjugate(this->ueem[s]);
|
||||
leemc[s] = conjugate(this->leem[s]);
|
||||
}
|
||||
deec[Ls] = conjugate(this->dee[Ls]);
|
||||
|
||||
this->MooeeInvCalls++;
|
||||
this->MooeeInvTime -= usecond();
|
||||
auto nloop = grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp1,tmp2;
|
||||
auto ss=sss*Ls;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// Apply (U^{\prime})^{-dagger}
|
||||
coalescedWrite(chi[ss], psi(ss));
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5m(tmp1, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s) - ueec[s-1]*tmp1);
|
||||
// Apply (U^{\prime})^{-dagger} U_m^{-\dagger}
|
||||
res = psi(ss);
|
||||
spProj5p(tmp,res);
|
||||
acc = conjugate(pueem[0])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= conjugate(puee[s-1])*tmp;
|
||||
spProj5p(tmp,res);
|
||||
acc += conjugate(pueem[s])*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1) - conjugate(puee[Ls-2])*tmp - acc;
|
||||
|
||||
// U_m^{-\dagger}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5p(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - ueemc[s]*tmp1);
|
||||
}
|
||||
|
||||
// L_m^{-\dagger} D^{-dagger}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5m(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s] ,(1.0/deec[s])*chi(ss+s) - (leemc[s]/deec[Ls-1])*tmp1);
|
||||
}
|
||||
spProj5p(tmp2, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/deec[Ls-1])*tmp1 + (1.0/deec[Ls])*tmp2);
|
||||
|
||||
// Apply L^{-dagger}
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
spProj5p(tmp1, chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s],chi(ss+s) - leec[s]*tmp1);
|
||||
// Apply L_m^{-\dagger} D^{-dagger} L^{-dagger}
|
||||
acc = conjugate(1.0/pdee[Ls-1])*res;
|
||||
tmp = conjugate(1.0/pdee[Ls ])*res;
|
||||
spProj5m(acc,acc);
|
||||
spProj5p(tmp,tmp);
|
||||
coalescedWrite(chi[ss+Ls-1], acc + tmp);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = conjugate(1.0/pdee[s])*chi(ss+s) - conjugate(plee[s])*tmp - conjugate(pleem[s])*acc;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
|
@ -11,6 +11,7 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
Author: Peter Boyle <peterboyle@Peters-MacBook-Pro-2.local>
|
||||
Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
Author: David Murphy <dmurphy@phys.columbia.edu>
|
||||
Author: Gianluca Filaci <g.filaci@ed.ac.uk>
|
||||
|
||||
This program is free software; you can redistribute it and/or modify
|
||||
it under the terms of the GNU General Public License as published by
|
||||
@ -49,6 +50,10 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField
|
||||
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
this->M5Dtime -= usecond();
|
||||
@ -64,7 +69,7 @@ void MobiusEOFAFermion<Impl>::M5D(const FermionField &psi_i, const FermionField
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5m(tmp1, psi(idx_u));
|
||||
spProj5p(tmp2, psi(idx_l));
|
||||
coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
|
||||
@ -88,6 +93,11 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const Fermion
|
||||
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
auto pshift_coeffs = &shift_coeffs[0];
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
this->M5Dtime -= usecond();
|
||||
@ -108,7 +118,7 @@ void MobiusEOFAFermion<Impl>::M5D_shift(const FermionField &psi_i, const Fermion
|
||||
if(pm == 1){ spProj5p(tmp, psi(ss+shift_s)); }
|
||||
else { spProj5m(tmp, psi(ss+shift_s)); }
|
||||
|
||||
coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 +lower[s]*tmp2 + shift_coeffs[s]*tmp);
|
||||
coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 +plower[s]*tmp2 + pshift_coeffs[s]*tmp);
|
||||
}
|
||||
});
|
||||
|
||||
@ -128,6 +138,10 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionFie
|
||||
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
this->M5Dtime -= usecond();
|
||||
@ -144,7 +158,7 @@ void MobiusEOFAFermion<Impl>::M5Ddag(const FermionField &psi_i, const FermionFie
|
||||
uint64_t idx_l = ss+((s+Ls-1)%Ls);
|
||||
spProj5p(tmp1, psi(idx_u));
|
||||
spProj5m(tmp2, psi(idx_l));
|
||||
coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
}
|
||||
});
|
||||
|
||||
@ -166,6 +180,11 @@ void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi_i, const Ferm
|
||||
|
||||
assert(phi.Checkerboard() == psi.Checkerboard());
|
||||
|
||||
auto pdiag = &diag[0];
|
||||
auto pupper = &upper[0];
|
||||
auto plower = &lower[0];
|
||||
auto pshift_coeffs = &shift_coeffs[0];
|
||||
|
||||
// Flops = 6.0*(Nc*Ns) *Ls*vol
|
||||
this->M5Dcalls++;
|
||||
this->M5Dtime -= usecond();
|
||||
@ -189,12 +208,12 @@ void MobiusEOFAFermion<Impl>::M5Ddag_shift(const FermionField &psi_i, const Ferm
|
||||
spProj5p(tmp1, psi(idx_u));
|
||||
spProj5m(tmp2, psi(idx_l));
|
||||
|
||||
if(s==(Ls-1)) coalescedWrite(chi[ss+s], chi(ss+s)+ diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
else coalescedWrite(chi[ss+s], diag[s]*phi(ss+s) + upper[s]*tmp1 + lower[s]*tmp2);
|
||||
if(s==(Ls-1)) coalescedWrite(chi[ss+s], chi(ss+s)+ pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
else coalescedWrite(chi[ss+s], pdiag[s]*phi(ss+s) + pupper[s]*tmp1 + plower[s]*tmp2);
|
||||
if(pm == 1){ spProj5p(tmp, psi(ss+s)); }
|
||||
else { spProj5m(tmp, psi(ss+s)); }
|
||||
|
||||
coalescedWrite(chi[ss+shift_s],chi(ss+shift_s)+shift_coeffs[s]*tmp);
|
||||
coalescedWrite(chi[ss+shift_s],chi(ss+shift_s)+pshift_coeffs[s]*tmp);
|
||||
}
|
||||
});
|
||||
|
||||
@ -223,36 +242,38 @@ void MobiusEOFAFermion<Impl>::MooeeInv(const FermionField &psi_i, FermionField &
|
||||
|
||||
int nloop = grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
|
||||
uint64_t ss = sss*Ls;
|
||||
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// Apply (L^{\prime})^{-1}
|
||||
coalescedWrite(chi[ss], psi(ss)); // chi[0]=psi[0]
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5p(tmp, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s) - plee[s-1]*tmp);
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (L^{\prime})^{-1} L_m^{-1}
|
||||
res = psi(ss);
|
||||
spProj5m(tmp,res);
|
||||
acc = pleem[0]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= plee[s-1]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
acc += pleem[s]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1) - plee[Ls-2]*tmp - acc;
|
||||
|
||||
// L_m^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
|
||||
spProj5m(tmp, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pleem[s]*tmp);
|
||||
}
|
||||
|
||||
// U_m^{-1} D^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[s] + 1/d chi[s]
|
||||
spProj5p(tmp, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s) - (pueem[s]/pdee[Ls-1])*tmp);
|
||||
}
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
|
||||
// Apply U^{-1}
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
spProj5m(tmp, chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - puee[s]*tmp);
|
||||
// Apply U_m^{-1} D^{-1} U^{-1}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
coalescedWrite(chi[ss+Ls-1],res);
|
||||
spProj5p(acc,res);
|
||||
spProj5m(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = (1.0/pdee[s])*chi(ss+s) - puee[s]*tmp - pueem[s]*acc;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
@ -281,44 +302,44 @@ void MobiusEOFAFermion<Impl>::MooeeInv_shift(const FermionField &psi_i, FermionF
|
||||
|
||||
int nloop = grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
|
||||
uint64_t ss = sss*Ls;
|
||||
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp1,tmp2,tmp2_spProj;
|
||||
spinor tmp, acc, res, tmp_spProj;
|
||||
|
||||
// Apply (L^{\prime})^{-1} and accumulate MooeeInv_shift_lc[j]*psi[j] in tmp2
|
||||
coalescedWrite(chi[ss], psi(ss)); // chi[0]=psi[0]
|
||||
tmp2 = pMooeeInv_shift_lc[0]*psi(ss);
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5p(tmp1, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s) - plee[s-1]*tmp1);
|
||||
tmp2 = tmp2 + pMooeeInv_shift_lc[s]*psi(ss+s);
|
||||
// Apply (L^{\prime})^{-1} L_m^{-1}
|
||||
res = psi(ss);
|
||||
spProj5m(tmp,res);
|
||||
acc = pleem[0]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
tmp_spProj = pMooeeInv_shift_lc[0]*res;
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
tmp_spProj += pMooeeInv_shift_lc[s]*res;
|
||||
res -= plee[s-1]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
acc += pleem[s]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
if(pm == 1){ spProj5p(tmp2_spProj, tmp2);}
|
||||
else { spProj5m(tmp2_spProj, tmp2); }
|
||||
res = psi(ss+Ls-1);
|
||||
|
||||
// L_m^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[ee] = 1 - sum[s<Ls-1] -leem[s]P_- chi
|
||||
spProj5m(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pleem[s]*tmp1);
|
||||
}
|
||||
tmp_spProj += pMooeeInv_shift_lc[Ls-1]*res;
|
||||
if(pm == 1){ spProj5p(tmp_spProj, tmp_spProj);}
|
||||
else { spProj5m(tmp_spProj, tmp_spProj); }
|
||||
|
||||
// U_m^{-1} D^{-1}
|
||||
for(int s=0; s<Ls-1; s++){ // Chi[s] + 1/d chi[s]
|
||||
spProj5p(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s) - (pueem[s]/pdee[Ls-1])*tmp1);
|
||||
}
|
||||
// chi[ss+Ls-1] = (1.0/pdee[Ls-1])*chi[ss+Ls-1] + MooeeInv_shift_norm[Ls-1]*tmp2_spProj;
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
spProj5m(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) + pMooeeInv_shift_norm[Ls-1]*tmp2_spProj);
|
||||
res = res - plee[Ls-2]*tmp - acc;
|
||||
|
||||
// Apply U^{-1} and add shift term
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
coalescedWrite(chi[ss+s] , chi(ss+s) - puee[s]*tmp1);
|
||||
spProj5m(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) + pMooeeInv_shift_norm[s]*tmp2_spProj);
|
||||
// Apply U_m^{-1} D^{-1} U^{-1}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
spProj5p(acc,res);
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+Ls-1], res + pMooeeInv_shift_norm[Ls-1]*tmp_spProj);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = (1.0/pdee[s])*chi(ss+s) - puee[s]*tmp - pueem[s]*acc;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s], res + pMooeeInv_shift_norm[s]*tmp_spProj);
|
||||
}
|
||||
});
|
||||
|
||||
@ -347,39 +368,40 @@ void MobiusEOFAFermion<Impl>::MooeeInvDag(const FermionField &psi_i, FermionFiel
|
||||
|
||||
int nloop = grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
|
||||
uint64_t ss = sss*Ls;
|
||||
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp;
|
||||
spinor tmp, acc, res;
|
||||
|
||||
// Apply (U^{\prime})^{-dag}
|
||||
coalescedWrite(chi[ss], psi(ss));
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5m(tmp, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s], psi(ss+s) - puee[s-1]*tmp);
|
||||
// X = Nc*Ns
|
||||
// flops = 2X + (Ls-2)(4X + 4X) + 6X + 1 + 2X + (Ls-1)(10X + 1) = -16X + Ls(1+18X) = -192 + 217*Ls flops
|
||||
// Apply (U^{\prime})^{-dagger} U_m^{-\dagger}
|
||||
res = psi(ss);
|
||||
spProj5p(tmp,res);
|
||||
acc = pueem[0]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
res -= puee[s-1]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
acc += pueem[s]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1) - puee[Ls-2]*tmp - acc;
|
||||
|
||||
// U_m^{-\dag}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5p(tmp, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pueem[s]*tmp);
|
||||
}
|
||||
|
||||
// L_m^{-\dag} D^{-dag}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5m(tmp, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s) - (pleem[s]/pdee[Ls-1])*tmp);
|
||||
}
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
|
||||
// Apply L^{-dag}
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
spProj5p(tmp, chi(ss+s+1));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - plee[s]*tmp);
|
||||
// Apply L_m^{-\dagger} D^{-dagger} L^{-dagger}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
coalescedWrite(chi[ss+Ls-1],res);
|
||||
spProj5m(acc,res);
|
||||
spProj5p(tmp,res);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = (1.0/pdee[s])*chi(ss+s) - plee[s]*tmp - pleem[s]*acc;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
});
|
||||
|
||||
this->MooeeInvTime += usecond();
|
||||
}
|
||||
|
||||
@ -406,44 +428,44 @@ void MobiusEOFAFermion<Impl>::MooeeInvDag_shift(const FermionField &psi_i, Fermi
|
||||
|
||||
int nloop = grid->oSites()/Ls;
|
||||
accelerator_for(sss,nloop,Simd::Nsimd(),{
|
||||
|
||||
uint64_t ss = sss*Ls;
|
||||
|
||||
uint64_t ss=sss*Ls;
|
||||
typedef decltype(coalescedRead(psi[0])) spinor;
|
||||
spinor tmp1,tmp2,tmp2_spProj;
|
||||
spinor tmp, acc, res, tmp_spProj;
|
||||
|
||||
// Apply (U^{\prime})^{-dag} and accumulate MooeeInvDag_shift_lc[j]*psi[j] in tmp2
|
||||
coalescedWrite(chi[ss], psi(ss));
|
||||
tmp2 = pMooeeInvDag_shift_lc[0]*psi(ss);
|
||||
for(int s=1; s<Ls; s++){
|
||||
spProj5m(tmp1, chi(ss+s-1));
|
||||
coalescedWrite(chi[ss+s],psi(ss+s) - puee[s-1]*tmp1);
|
||||
tmp2 = tmp2 + pMooeeInvDag_shift_lc[s]*psi(ss+s);
|
||||
// Apply (U^{\prime})^{-dagger} U_m^{-\dagger}
|
||||
res = psi(ss);
|
||||
spProj5p(tmp,res);
|
||||
acc = pueem[0]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss],res);
|
||||
tmp_spProj = pMooeeInvDag_shift_lc[0]*res;
|
||||
|
||||
for(int s=1;s<Ls-1;s++){
|
||||
res = psi(ss+s);
|
||||
tmp_spProj += pMooeeInvDag_shift_lc[s]*res;
|
||||
res -= puee[s-1]*tmp;
|
||||
spProj5p(tmp,res);
|
||||
acc += pueem[s]*tmp;
|
||||
spProj5m(tmp,res);
|
||||
coalescedWrite(chi[ss+s],res);
|
||||
}
|
||||
res = psi(ss+Ls-1);
|
||||
|
||||
if(pm == 1){ spProj5p(tmp2_spProj, tmp2);}
|
||||
else { spProj5m(tmp2_spProj, tmp2);}
|
||||
tmp_spProj += pMooeeInvDag_shift_lc[Ls-1]*res;
|
||||
if(pm == 1){ spProj5p(tmp_spProj, tmp_spProj); }
|
||||
else { spProj5m(tmp_spProj, tmp_spProj); }
|
||||
|
||||
// U_m^{-\dag}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5p(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) - pueem[s]*tmp1);
|
||||
}
|
||||
res = res - puee[Ls-2]*tmp - acc;
|
||||
|
||||
// L_m^{-\dag} D^{-dag}
|
||||
for(int s=0; s<Ls-1; s++){
|
||||
spProj5m(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+s], (1.0/pdee[s])*chi(ss+s) - (pleem[s]/pdee[Ls-1])*tmp1);
|
||||
}
|
||||
coalescedWrite(chi[ss+Ls-1], (1.0/pdee[Ls-1])*chi(ss+Ls-1));
|
||||
spProj5p(tmp1, chi(ss+Ls-1));
|
||||
coalescedWrite(chi[ss+Ls-1], chi(ss+Ls-1) + pMooeeInvDag_shift_norm[Ls-1]*tmp2_spProj);
|
||||
|
||||
// Apply L^{-dag}
|
||||
for(int s=Ls-2; s>=0; s--){
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) - plee[s]*tmp1);
|
||||
spProj5p(tmp1, chi(ss+s));
|
||||
coalescedWrite(chi[ss+s], chi(ss+s) + pMooeeInvDag_shift_norm[s]*tmp2_spProj);
|
||||
// Apply L_m^{-\dagger} D^{-dagger} L^{-dagger}
|
||||
res = (1.0/pdee[Ls-1])*res;
|
||||
spProj5m(acc,res);
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+Ls-1], res + pMooeeInvDag_shift_norm[Ls-1]*tmp_spProj);
|
||||
for (int s=Ls-2;s>=0;s--){
|
||||
res = (1.0/pdee[s])*chi(ss+s) - plee[s]*tmp - pleem[s]*acc;
|
||||
spProj5p(tmp,res);
|
||||
coalescedWrite(chi[ss+s], res + pMooeeInvDag_shift_norm[s]*tmp_spProj);
|
||||
}
|
||||
});
|
||||
|
||||
|
@ -97,7 +97,6 @@ protected:
|
||||
////////////////////////////////////
|
||||
// Classes for the user
|
||||
////////////////////////////////////
|
||||
// Note: the space time grid should be out of the QCD namespace
|
||||
template <class vector_type>
|
||||
class GridFourDimModule : public GridModule
|
||||
{
|
||||
|
@ -27,8 +27,7 @@ directory
|
||||
*************************************************************************************/
|
||||
#pragma once
|
||||
|
||||
namespace Grid {
|
||||
namespace QCD {
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
template <class Gimpl> class CovariantSmearing : public Gimpl
|
||||
{
|
||||
@ -84,4 +83,5 @@ public:
|
||||
}
|
||||
}
|
||||
};
|
||||
}}
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
@ -201,7 +201,6 @@ void G5R5(Lattice<vobj> &z,const Lattice<vobj> &x)
|
||||
});
|
||||
}
|
||||
|
||||
// I explicitly need these outside the QCD namespace
|
||||
template<typename vobj>
|
||||
void G5C(Lattice<vobj> &z, const Lattice<vobj> &x)
|
||||
{
|
||||
|
@ -31,7 +31,6 @@ directory
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
int threads = GridThread::GetThreads();
|
||||
|
@ -34,8 +34,7 @@ directory
|
||||
#define MIXED_PRECISION
|
||||
#endif
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
/*
|
||||
* Need a plan for gauge field update for mixed precision in HMC (2x speed up)
|
||||
@ -163,11 +162,11 @@ namespace Grid{
|
||||
MPCG(src,psi);
|
||||
}
|
||||
};
|
||||
}};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
int threads = GridThread::GetThreads();
|
||||
|
@ -34,8 +34,7 @@ directory
|
||||
#define MIXED_PRECISION
|
||||
#endif
|
||||
|
||||
namespace Grid{
|
||||
namespace QCD{
|
||||
NAMESPACE_BEGIN(Grid);
|
||||
|
||||
/*
|
||||
* Need a plan for gauge field update for mixed precision in HMC (2x speed up)
|
||||
@ -146,11 +145,12 @@ namespace Grid{
|
||||
MPCG(src,psi);
|
||||
}
|
||||
};
|
||||
}};
|
||||
|
||||
NAMESPACE_END(Grid);
|
||||
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
int threads = GridThread::GetThreads();
|
||||
|
@ -31,7 +31,6 @@ directory
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
int threads = GridThread::GetThreads();
|
||||
|
@ -1,8 +1,8 @@
|
||||
#include <Hadrons/Modules/MAction/DWF.hpp>
|
||||
#include <Hadrons/Modules/MAction/MobiusDWF.hpp>
|
||||
#include <Hadrons/Modules/MAction/ScaledDWF.hpp>
|
||||
#include <Hadrons/Modules/MAction/Wilson.hpp>
|
||||
#include <Hadrons/Modules/MAction/WilsonClover.hpp>
|
||||
#include <Hadrons/Modules/MAction/Wilson.hpp>
|
||||
#include <Hadrons/Modules/MAction/ZMobiusDWF.hpp>
|
||||
#include <Hadrons/Modules/MContraction/A2AAslashField.hpp>
|
||||
#include <Hadrons/Modules/MContraction/A2AFourQuarkContraction.hpp>
|
||||
@ -33,8 +33,8 @@
|
||||
#include <Hadrons/Modules/MGauge/Random.hpp>
|
||||
#include <Hadrons/Modules/MGauge/StochEm.hpp>
|
||||
#include <Hadrons/Modules/MGauge/StoutSmearing.hpp>
|
||||
#include <Hadrons/Modules/MGauge/Unit.hpp>
|
||||
#include <Hadrons/Modules/MGauge/UnitEm.hpp>
|
||||
#include <Hadrons/Modules/MGauge/Unit.hpp>
|
||||
#include <Hadrons/Modules/MIO/LoadA2AMatrixDiskVector.hpp>
|
||||
#include <Hadrons/Modules/MIO/LoadA2AVectors.hpp>
|
||||
#include <Hadrons/Modules/MIO/LoadBinary.hpp>
|
||||
@ -44,12 +44,12 @@
|
||||
#include <Hadrons/Modules/MIO/LoadEigenPack.hpp>
|
||||
#include <Hadrons/Modules/MIO/LoadNersc.hpp>
|
||||
#include <Hadrons/Modules/MIO/LoadPerambulator.hpp>
|
||||
#include <Hadrons/Modules/MNPR/Amputate.hpp>
|
||||
#include <Hadrons/Modules/MNPR/Bilinear.hpp>
|
||||
#include <Hadrons/Modules/MNPR/FourQuark.hpp>
|
||||
#include <Hadrons/Modules/MNoise/FullVolumeSpinColorDiagonal.hpp>
|
||||
#include <Hadrons/Modules/MNoise/SparseSpinColorDiagonal.hpp>
|
||||
#include <Hadrons/Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp>
|
||||
#include <Hadrons/Modules/MNPR/Amputate.hpp>
|
||||
#include <Hadrons/Modules/MNPR/Bilinear.hpp>
|
||||
#include <Hadrons/Modules/MNPR/FourQuark.hpp>
|
||||
#include <Hadrons/Modules/MScalar/ChargedProp.hpp>
|
||||
#include <Hadrons/Modules/MScalar/FreeProp.hpp>
|
||||
#include <Hadrons/Modules/MScalar/Scalar.hpp>
|
||||
@ -57,10 +57,10 @@
|
||||
#include <Hadrons/Modules/MScalarSUN/EMT.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/Grad.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/StochFreeField.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TrKinetic.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TrMag.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TrPhi.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TransProj.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TwoPoint.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/TwoPointNPR.hpp>
|
||||
#include <Hadrons/Modules/MScalarSUN/Utils.hpp>
|
||||
@ -74,6 +74,7 @@
|
||||
#include <Hadrons/Modules/MSolver/RBPrecCG.hpp>
|
||||
#include <Hadrons/Modules/MSource/Convolution.hpp>
|
||||
#include <Hadrons/Modules/MSource/Gauss.hpp>
|
||||
#include <Hadrons/Modules/MSource/JacobiSmear.hpp>
|
||||
#include <Hadrons/Modules/MSource/Momentum.hpp>
|
||||
#include <Hadrons/Modules/MSource/MomentumPhase.hpp>
|
||||
#include <Hadrons/Modules/MSource/Point.hpp>
|
||||
|
7
Hadrons/Modules/MSource/JacobiSmear.cc
Normal file
7
Hadrons/Modules/MSource/JacobiSmear.cc
Normal file
@ -0,0 +1,7 @@
|
||||
#include <Hadrons/Modules/MSource/JacobiSmear.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Hadrons;
|
||||
using namespace MSource;
|
||||
|
||||
template class Grid::Hadrons::MSource::TJacobiSmear<FIMPL>;
|
105
Hadrons/Modules/MSource/JacobiSmear.hpp
Normal file
105
Hadrons/Modules/MSource/JacobiSmear.hpp
Normal file
@ -0,0 +1,105 @@
|
||||
#ifndef Hadrons_MSource_JacobiSmear_hpp_
|
||||
#define Hadrons_MSource_JacobiSmear_hpp_
|
||||
|
||||
#include <Hadrons/Global.hpp>
|
||||
#include <Hadrons/Module.hpp>
|
||||
#include <Hadrons/ModuleFactory.hpp>
|
||||
|
||||
BEGIN_HADRONS_NAMESPACE
|
||||
|
||||
/******************************************************************************
|
||||
* JacobiSmear *
|
||||
******************************************************************************/
|
||||
BEGIN_MODULE_NAMESPACE(MSource)
|
||||
|
||||
class JacobiSmearPar: Serializable
|
||||
{
|
||||
public:
|
||||
GRID_SERIALIZABLE_CLASS_MEMBERS(JacobiSmearPar,
|
||||
std::string, gauge,
|
||||
double, width,
|
||||
int, iterations,
|
||||
int, orthog,
|
||||
std::string, source);
|
||||
};
|
||||
|
||||
template <typename FImpl>
|
||||
class TJacobiSmear: public Module<JacobiSmearPar>
|
||||
{
|
||||
public:
|
||||
FERM_TYPE_ALIASES(FImpl,);
|
||||
typedef typename FImpl::GaugeLinkField GaugeMat;
|
||||
public:
|
||||
// constructor
|
||||
TJacobiSmear(const std::string name);
|
||||
// destructor
|
||||
virtual ~TJacobiSmear(void) {};
|
||||
// dependency relation
|
||||
virtual std::vector<std::string> getInput(void);
|
||||
virtual std::vector<std::string> getOutput(void);
|
||||
// setup
|
||||
virtual void setup(void);
|
||||
// execution
|
||||
virtual void execute(void);
|
||||
};
|
||||
|
||||
MODULE_REGISTER_TMP(JacobiSmear, TJacobiSmear<FIMPL>, MSource);
|
||||
|
||||
/******************************************************************************
|
||||
* TJacobiSmear implementation *
|
||||
******************************************************************************/
|
||||
// constructor /////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
TJacobiSmear<FImpl>::TJacobiSmear(const std::string name)
|
||||
: Module<JacobiSmearPar>(name)
|
||||
{}
|
||||
|
||||
// dependencies/products ///////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TJacobiSmear<FImpl>::getInput(void)
|
||||
{
|
||||
std::vector<std::string> in = {par().source, par().gauge};
|
||||
|
||||
return in;
|
||||
}
|
||||
|
||||
template <typename FImpl>
|
||||
std::vector<std::string> TJacobiSmear<FImpl>::getOutput(void)
|
||||
{
|
||||
std::vector<std::string> out = {getName()};
|
||||
|
||||
return out;
|
||||
}
|
||||
|
||||
// setup ///////////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TJacobiSmear<FImpl>::setup(void)
|
||||
{
|
||||
envCreateLat(PropagatorField, getName());
|
||||
envTmp(std::vector<GaugeMat>, "Umu", 1, 4, envGetGrid(LatticeColourMatrix));
|
||||
}
|
||||
|
||||
// execution ///////////////////////////////////////////////////////////////////
|
||||
template <typename FImpl>
|
||||
void TJacobiSmear<FImpl>::execute(void)
|
||||
{
|
||||
auto &out = envGet(PropagatorField, getName());
|
||||
auto &src = envGet(PropagatorField, par().source);
|
||||
auto &U = envGet(GaugeField, par().gauge);
|
||||
envGetTmp(std::vector<GaugeMat>, Umu);
|
||||
for(int mu=0; mu<4; mu++)
|
||||
{
|
||||
Umu.at(mu)=peekLorentz(U,mu);
|
||||
}
|
||||
CovariantSmearing<FImpl> covsmear;
|
||||
out=src;
|
||||
startTimer("Jacobi iteration");
|
||||
covsmear.GaussianSmear(Umu, out, par().width, par().iterations, par().orthog);
|
||||
stopTimer("Jacobi iteration");
|
||||
}
|
||||
|
||||
END_MODULE_NAMESPACE
|
||||
|
||||
END_HADRONS_NAMESPACE
|
||||
|
||||
#endif // Hadrons_MSource_JacobiSmear_hpp_
|
@ -29,7 +29,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
|
||||
#include <Hadrons/Application.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace QCD;
|
||||
using namespace Hadrons;
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
|
@ -44,22 +44,22 @@ modules_cc =\
|
||||
Modules/MIO/LoadEigenPack.cc \
|
||||
Modules/MIO/LoadNersc.cc \
|
||||
Modules/MIO/LoadPerambulator.cc \
|
||||
Modules/MNPR/Amputate.cc \
|
||||
Modules/MNPR/Bilinear.cc \
|
||||
Modules/MNPR/FourQuark.cc \
|
||||
Modules/MNoise/FullVolumeSpinColorDiagonal.cc \
|
||||
Modules/MNoise/SparseSpinColorDiagonal.cc \
|
||||
Modules/MNoise/TimeDilutedSpinColorDiagonal.cc \
|
||||
Modules/MNPR/Amputate.cc \
|
||||
Modules/MNPR/Bilinear.cc \
|
||||
Modules/MNPR/FourQuark.cc \
|
||||
Modules/MScalar/ChargedProp.cc \
|
||||
Modules/MScalar/FreeProp.cc \
|
||||
Modules/MScalarSUN/Div.cc \
|
||||
Modules/MScalarSUN/EMT.cc \
|
||||
Modules/MScalarSUN/Grad.cc \
|
||||
Modules/MScalarSUN/StochFreeField.cc \
|
||||
Modules/MScalarSUN/TransProj.cc \
|
||||
Modules/MScalarSUN/TrKinetic.cc \
|
||||
Modules/MScalarSUN/TrMag.cc \
|
||||
Modules/MScalarSUN/TrPhi.cc \
|
||||
Modules/MScalarSUN/TransProj.cc \
|
||||
Modules/MScalarSUN/TwoPoint.cc \
|
||||
Modules/MScalarSUN/TwoPointNPR.cc \
|
||||
Modules/MSink/Point.cc \
|
||||
@ -71,6 +71,7 @@ modules_cc =\
|
||||
Modules/MSolver/RBPrecCG.cc \
|
||||
Modules/MSource/Convolution.cc \
|
||||
Modules/MSource/Gauss.cc \
|
||||
Modules/MSource/JacobiSmear.cc \
|
||||
Modules/MSource/Momentum.cc \
|
||||
Modules/MSource/MomentumPhase.cc \
|
||||
Modules/MSource/Point.cc \
|
||||
@ -86,8 +87,8 @@ modules_hpp =\
|
||||
Modules/MAction/DWF.hpp \
|
||||
Modules/MAction/MobiusDWF.hpp \
|
||||
Modules/MAction/ScaledDWF.hpp \
|
||||
Modules/MAction/Wilson.hpp \
|
||||
Modules/MAction/WilsonClover.hpp \
|
||||
Modules/MAction/Wilson.hpp \
|
||||
Modules/MAction/ZMobiusDWF.hpp \
|
||||
Modules/MContraction/A2AAslashField.hpp \
|
||||
Modules/MContraction/A2AFourQuarkContraction.hpp \
|
||||
@ -118,8 +119,8 @@ modules_hpp =\
|
||||
Modules/MGauge/Random.hpp \
|
||||
Modules/MGauge/StochEm.hpp \
|
||||
Modules/MGauge/StoutSmearing.hpp \
|
||||
Modules/MGauge/Unit.hpp \
|
||||
Modules/MGauge/UnitEm.hpp \
|
||||
Modules/MGauge/Unit.hpp \
|
||||
Modules/MIO/LoadA2AMatrixDiskVector.hpp \
|
||||
Modules/MIO/LoadA2AVectors.hpp \
|
||||
Modules/MIO/LoadBinary.hpp \
|
||||
@ -129,12 +130,12 @@ modules_hpp =\
|
||||
Modules/MIO/LoadEigenPack.hpp \
|
||||
Modules/MIO/LoadNersc.hpp \
|
||||
Modules/MIO/LoadPerambulator.hpp \
|
||||
Modules/MNPR/Amputate.hpp \
|
||||
Modules/MNPR/Bilinear.hpp \
|
||||
Modules/MNPR/FourQuark.hpp \
|
||||
Modules/MNoise/FullVolumeSpinColorDiagonal.hpp \
|
||||
Modules/MNoise/SparseSpinColorDiagonal.hpp \
|
||||
Modules/MNoise/TimeDilutedSpinColorDiagonal.hpp \
|
||||
Modules/MNPR/Amputate.hpp \
|
||||
Modules/MNPR/Bilinear.hpp \
|
||||
Modules/MNPR/FourQuark.hpp \
|
||||
Modules/MScalar/ChargedProp.hpp \
|
||||
Modules/MScalar/FreeProp.hpp \
|
||||
Modules/MScalar/Scalar.hpp \
|
||||
@ -142,10 +143,10 @@ modules_hpp =\
|
||||
Modules/MScalarSUN/EMT.hpp \
|
||||
Modules/MScalarSUN/Grad.hpp \
|
||||
Modules/MScalarSUN/StochFreeField.hpp \
|
||||
Modules/MScalarSUN/TransProj.hpp \
|
||||
Modules/MScalarSUN/TrKinetic.hpp \
|
||||
Modules/MScalarSUN/TrMag.hpp \
|
||||
Modules/MScalarSUN/TrPhi.hpp \
|
||||
Modules/MScalarSUN/TransProj.hpp \
|
||||
Modules/MScalarSUN/TwoPoint.hpp \
|
||||
Modules/MScalarSUN/TwoPointNPR.hpp \
|
||||
Modules/MScalarSUN/Utils.hpp \
|
||||
@ -159,6 +160,7 @@ modules_hpp =\
|
||||
Modules/MSolver/RBPrecCG.hpp \
|
||||
Modules/MSource/Convolution.hpp \
|
||||
Modules/MSource/Gauss.hpp \
|
||||
Modules/MSource/JacobiSmear.hpp \
|
||||
Modules/MSource/Momentum.hpp \
|
||||
Modules/MSource/MomentumPhase.hpp \
|
||||
Modules/MSource/Point.hpp \
|
||||
|
13
bootstrap.sh
13
bootstrap.sh
@ -1,10 +1,21 @@
|
||||
#!/usr/bin/env bash
|
||||
set -e
|
||||
|
||||
EIGEN_URL='https://gitlab.com/libeigen/eigen/-/archive/3.3.7/eigen-3.3.7.tar.bz2'
|
||||
EIGEN_SHA256SUM='685adf14bd8e9c015b78097c1dc22f2f01343756f196acdc76a678e1ae352e11'
|
||||
|
||||
|
||||
echo "-- deploying Eigen source..."
|
||||
ARC=`basename ${EIGEN_URL}`
|
||||
wget ${EIGEN_URL} --no-check-certificate && ./scripts/update_eigen.sh ${ARC} && rm ${ARC}
|
||||
wget ${EIGEN_URL} --no-check-certificate
|
||||
if command -v sha256sum; then
|
||||
echo "$EIGEN_SHA256SUM $(basename "$EIGEN_URL")" \
|
||||
| sha256sum --check || exit 1
|
||||
else
|
||||
echo "WARNING: could not verify checksum, please install sha256sum" >&2
|
||||
fi
|
||||
./scripts/update_eigen.sh ${ARC}
|
||||
rm ${ARC}
|
||||
# patch for non-portable includes in Eigen 3.3.5
|
||||
# apparently already fixed in Eigen HEAD so it should not be
|
||||
# a problem in the future (A.P.)
|
||||
|
@ -281,6 +281,9 @@ case ${CXX} in
|
||||
CXX="nvcc -x cu "
|
||||
CXXLD="nvcc -link"
|
||||
CXXFLAGS="$CXXFLAGS -Xcompiler -fno-strict-aliasing -Xcompiler -Wno-unusable-partial-specialization --expt-extended-lambda --expt-relaxed-constexpr"
|
||||
if test $ac_openmp = yes; then
|
||||
CXXFLAGS="$CXXFLAGS -Xcompiler -fopenmp"
|
||||
fi
|
||||
;;
|
||||
*)
|
||||
CXXLD=${CXX}
|
||||
|
@ -28,7 +28,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
|
||||
#include <Hadrons/DiskVector.hpp>
|
||||
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
using namespace Grid::Hadrons;
|
||||
|
||||
GRID_SERIALIZABLE_ENUM(Enum, undef, red, 1, blue, 2, green, 3);
|
||||
|
@ -31,7 +31,6 @@ directory
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid_init(&argc, &argv);
|
||||
int threads = GridThread::GetThreads();
|
||||
|
@ -27,15 +27,15 @@ Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
|
||||
/* END LEGAL */
|
||||
#include <Grid/Grid.h>
|
||||
|
||||
double calc_grid_p (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_chroma_p (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_grid_r (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_grid_IW (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_grid_r_dir (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_chroma_r (Grid::QCD::LatticeGaugeField & lat);
|
||||
double calc_chroma_IW (Grid::QCD::LatticeGaugeField & lat);
|
||||
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu);
|
||||
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu);
|
||||
double calc_grid_p (Grid::LatticeGaugeField & lat);
|
||||
double calc_chroma_p (Grid::LatticeGaugeField & lat);
|
||||
double calc_grid_r (Grid::LatticeGaugeField & lat);
|
||||
double calc_grid_IW (Grid::LatticeGaugeField & lat);
|
||||
double calc_grid_r_dir (Grid::LatticeGaugeField & lat);
|
||||
double calc_chroma_r (Grid::LatticeGaugeField & lat);
|
||||
double calc_chroma_IW (Grid::LatticeGaugeField & lat);
|
||||
void check_grid_r_staple(Grid::LatticeGaugeField & Umu);
|
||||
void check_grid_p_staple(Grid::LatticeGaugeField & Umu);
|
||||
|
||||
const double beta=2.6;
|
||||
const double c1=-0.331;
|
||||
@ -53,10 +53,10 @@ public:
|
||||
|
||||
typedef multi1d<LatticeColorMatrix> U;
|
||||
|
||||
static void ImportGauge(Grid::QCD::LatticeGaugeField & gr,
|
||||
static void ImportGauge(Grid::LatticeGaugeField & gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -127,8 +127,8 @@ int main (int argc,char **argv )
|
||||
* Setup Grid
|
||||
*********************************************************/
|
||||
Grid::Grid_init(&argc,&argv);
|
||||
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridDefaultMpi());
|
||||
|
||||
std::vector<int> gd = UGrid->GlobalDimensions();
|
||||
@ -138,7 +138,7 @@ int main (int argc,char **argv )
|
||||
QDP::Layout::setLattSize(nrow);
|
||||
QDP::Layout::create();
|
||||
|
||||
Grid::QCD::LatticeGaugeField lat(UGrid);
|
||||
Grid::LatticeGaugeField lat(UGrid);
|
||||
|
||||
double s_grid = calc_grid_p (lat);
|
||||
|
||||
@ -181,7 +181,7 @@ int main (int argc,char **argv )
|
||||
Chroma::finalize();
|
||||
}
|
||||
|
||||
double calc_chroma_IW(Grid::QCD::LatticeGaugeField & lat)
|
||||
double calc_chroma_IW(Grid::LatticeGaugeField & lat)
|
||||
{
|
||||
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
|
||||
|
||||
@ -203,7 +203,7 @@ double calc_chroma_IW(Grid::QCD::LatticeGaugeField & lat)
|
||||
|
||||
return s;
|
||||
}
|
||||
double calc_chroma_r(Grid::QCD::LatticeGaugeField & lat)
|
||||
double calc_chroma_r(Grid::LatticeGaugeField & lat)
|
||||
{
|
||||
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
|
||||
|
||||
@ -245,7 +245,7 @@ double calc_chroma_r(Grid::QCD::LatticeGaugeField & lat)
|
||||
// action = beta * Nd*Nd-1*vol*0.5 - beta * Nd*Nd-1*vol*0.5*plaq
|
||||
//
|
||||
// plaq == sumplaq * 2/(V*Nd*(Nd-1)*Nc)
|
||||
double calc_chroma_p(Grid::QCD::LatticeGaugeField & lat)
|
||||
double calc_chroma_p(Grid::LatticeGaugeField & lat)
|
||||
{
|
||||
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
|
||||
|
||||
@ -270,60 +270,60 @@ double calc_chroma_p(Grid::QCD::LatticeGaugeField & lat)
|
||||
|
||||
|
||||
|
||||
double calc_grid_p(Grid::QCD::LatticeGaugeField & Umu)
|
||||
double calc_grid_p(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
|
||||
Grid::SU3::HotConfiguration(RNG4,Umu);
|
||||
|
||||
Grid::QCD::LatticeColourMatrix tmp(UGrid);
|
||||
Grid::LatticeColourMatrix tmp(UGrid);
|
||||
tmp = Grid::zero;
|
||||
|
||||
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,2);
|
||||
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,3);
|
||||
Grid::PokeIndex<LorentzIndex>(Umu,tmp,2);
|
||||
Grid::PokeIndex<LorentzIndex>(Umu,tmp,3);
|
||||
|
||||
Grid::QCD::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
|
||||
Grid::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
|
||||
|
||||
return Wilson.S(Umu);
|
||||
}
|
||||
double calc_grid_r(Grid::QCD::LatticeGaugeField & Umu)
|
||||
double calc_grid_r(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
|
||||
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
|
||||
return Wilson.S(Umu);
|
||||
}
|
||||
double calc_grid_IW(Grid::QCD::LatticeGaugeField & Umu)
|
||||
double calc_grid_IW(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
|
||||
Grid::QCD::IwasakiGaugeActionR Iwasaki(beta);
|
||||
Grid::IwasakiGaugeActionR Iwasaki(beta);
|
||||
|
||||
return Iwasaki.S(Umu);
|
||||
}
|
||||
double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
double calc_grid_r_dir(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
|
||||
std::vector<Grid::QCD::LatticeColourMatrix> U(4,UGrid);
|
||||
std::vector<Grid::LatticeColourMatrix> U(4,UGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = Grid::PeekIndex<LorentzIndex>(Umu,mu);
|
||||
}
|
||||
|
||||
Grid::QCD::LatticeComplex rect(UGrid);
|
||||
Grid::QCD::TComplex trect;
|
||||
Grid::QCD::Complex crect;
|
||||
Grid::LatticeComplex rect(UGrid);
|
||||
Grid::TComplex trect;
|
||||
Grid::Complex crect;
|
||||
Grid::RealD rrect;
|
||||
Grid::RealD vol = UGrid->gSites();
|
||||
for(int mu=0;mu<Grid::QCD::Nd;mu++){
|
||||
for(int nu=0;nu<Grid::QCD::Nd;nu++){
|
||||
for(int mu=0;mu<Grid::Nd;mu++){
|
||||
for(int nu=0;nu<Grid::Nd;nu++){
|
||||
if ( mu!=nu ) {
|
||||
|
||||
Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
|
||||
Grid::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
|
||||
trect = Grid::sum(rect);
|
||||
crect = Grid::TensorRemove(trect);
|
||||
rrect = real(crect);
|
||||
@ -335,9 +335,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// Staple test
|
||||
Peter.Start();
|
||||
{
|
||||
Grid::QCD::LatticeColourMatrix Stap(UGrid);
|
||||
Grid::QCD::LatticeComplex SumTrStap(UGrid);
|
||||
Grid::QCD::LatticeComplex TrStap(UGrid);
|
||||
Grid::LatticeColourMatrix Stap(UGrid);
|
||||
Grid::LatticeComplex SumTrStap(UGrid);
|
||||
Grid::LatticeComplex TrStap(UGrid);
|
||||
|
||||
/*
|
||||
* Make staple for loops centered at coor of link ; this one is ok. // |
|
||||
@ -346,10 +346,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// __ ___
|
||||
// | __ |
|
||||
Stap =
|
||||
Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftForward (U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
@ -364,10 +364,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// __
|
||||
// |__ __ |
|
||||
|
||||
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
|
||||
Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
|
||||
@ -379,10 +379,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// __
|
||||
// |__ __ |
|
||||
|
||||
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
|
||||
Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
|
||||
@ -395,10 +395,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// __ ___
|
||||
// |__ |
|
||||
|
||||
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
|
||||
Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward (U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
|
||||
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
@ -418,12 +418,12 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
* Make staple for loops centered at coor of link ; this one is ok. // |
|
||||
*/
|
||||
// Stap =
|
||||
// Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
|
||||
// Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ;
|
||||
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
|
||||
// Grid::Cshift(Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
|
||||
// Grid::adj(Grid::PeriodicBC::CovShiftForward(U[nu],nu,Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ;
|
||||
Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftForward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftForward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
SumTrStap += TrStap;
|
||||
@ -440,10 +440,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// | |
|
||||
// --
|
||||
|
||||
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
|
||||
Stap = Grid::Cshift(Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[nu],nu,
|
||||
Grid::PeriodicBC::CovShiftBackward(U[mu],mu,
|
||||
Grid::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
|
||||
|
||||
TrStap = Grid::trace (U[mu]*Stap);
|
||||
trect = Grid::sum(TrStap);
|
||||
@ -459,20 +459,20 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
Peter.Stop();
|
||||
Azusa.Start();
|
||||
{
|
||||
Grid::QCD::LatticeComplex RectPlaq_d(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix ds_U(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix left_2(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix upper_l(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix upper_staple(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix down_l(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix down_staple(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix tmp(UGrid);
|
||||
Grid::LatticeComplex RectPlaq_d(UGrid);
|
||||
Grid::LatticeColourMatrix ds_U(UGrid);
|
||||
Grid::LatticeColourMatrix left_2(UGrid);
|
||||
Grid::LatticeColourMatrix upper_l(UGrid);
|
||||
Grid::LatticeColourMatrix upper_staple(UGrid);
|
||||
Grid::LatticeColourMatrix down_l(UGrid);
|
||||
Grid::LatticeColourMatrix down_staple(UGrid);
|
||||
Grid::LatticeColourMatrix tmp(UGrid);
|
||||
|
||||
// 2 (mu)x1(nu)
|
||||
left_2= Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
|
||||
left_2= Grid::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
|
||||
tmp=Grid::Cshift(U[nu],mu,2); // Unu(x+2mu)
|
||||
|
||||
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
|
||||
upper_l= Grid::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
|
||||
// __ __
|
||||
// = |
|
||||
|
||||
@ -546,9 +546,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// _
|
||||
// | |
|
||||
// | |
|
||||
Grid::QCD::LatticeColourMatrix up2= Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
|
||||
Grid::LatticeColourMatrix up2= Grid::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
|
||||
|
||||
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
|
||||
upper_l= Grid::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
|
||||
ds_U= upper_l*Grid::adj(up2);
|
||||
|
||||
RectPlaq_d = Grid::trace(U[mu]*ds_U);
|
||||
@ -569,7 +569,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
downer_l= |
|
||||
(x)<----V
|
||||
*/
|
||||
down_l= Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
|
||||
down_l= Grid::adj(Grid::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
|
||||
/*
|
||||
^ |
|
||||
down_staple = | V
|
||||
@ -601,23 +601,23 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
|
||||
|
||||
}
|
||||
}}
|
||||
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
|
||||
return Wilson.S(Umu);
|
||||
};
|
||||
|
||||
void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
|
||||
void check_grid_r_staple(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
|
||||
Grid::QCD::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
Grid::PlaqPlusRectangleActionR Wilson(0.0,c1); // Just take beta = 0.0
|
||||
|
||||
Grid::QCD::LatticeColourMatrix staple(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix link(UGrid);
|
||||
Grid::QCD::LatticeComplex Traced(UGrid);
|
||||
Grid::LatticeColourMatrix staple(UGrid);
|
||||
Grid::LatticeColourMatrix link(UGrid);
|
||||
Grid::LatticeComplex Traced(UGrid);
|
||||
Grid::Complex Rplaq(0.0);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
@ -630,12 +630,12 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// Vol as for each site
|
||||
Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0);
|
||||
|
||||
Grid::QCD::ColourWilsonLoops::RectStaple(staple,Umu,mu);
|
||||
Grid::ColourWilsonLoops::RectStaple(staple,Umu,mu);
|
||||
|
||||
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);
|
||||
link = Grid::PeekIndex<LorentzIndex>(Umu,mu);
|
||||
|
||||
Traced = Grid::trace( link*staple) * RectScale;
|
||||
Grid::QCD::TComplex Tp = Grid::sum(Traced);
|
||||
Grid::TComplex Tp = Grid::sum(Traced);
|
||||
Grid::Complex p = Grid::TensorRemove(Tp);
|
||||
|
||||
std::cout<< "Rect from RectStaple "<<p<<std::endl;
|
||||
@ -645,18 +645,18 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
|
||||
std::cout<< "Rect from RectStaple "<<Rplaq<<std::endl;
|
||||
}
|
||||
|
||||
void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
|
||||
void check_grid_p_staple(Grid::LatticeGaugeField & Umu)
|
||||
{
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
|
||||
Grid::QCD::PlaqPlusRectangleActionR Wilson(1.0,0.0); // Just take c1 = 0.0
|
||||
Grid::PlaqPlusRectangleActionR Wilson(1.0,0.0); // Just take c1 = 0.0
|
||||
|
||||
Grid::QCD::LatticeColourMatrix staple(UGrid);
|
||||
Grid::QCD::LatticeColourMatrix link(UGrid);
|
||||
Grid::QCD::LatticeComplex Traced(UGrid);
|
||||
Grid::LatticeColourMatrix staple(UGrid);
|
||||
Grid::LatticeColourMatrix link(UGrid);
|
||||
Grid::LatticeComplex Traced(UGrid);
|
||||
Grid::Complex plaq(0.0);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
@ -669,12 +669,12 @@ void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
|
||||
// Vol as for each site
|
||||
Grid::RealD Scale(1.0/vol/12.0/2.0/3.0);
|
||||
|
||||
Grid::QCD::ColourWilsonLoops::Staple(staple,Umu,mu);
|
||||
Grid::ColourWilsonLoops::Staple(staple,Umu,mu);
|
||||
|
||||
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);
|
||||
link = Grid::PeekIndex<LorentzIndex>(Umu,mu);
|
||||
|
||||
Traced = Grid::trace( link*staple) * Scale;
|
||||
Grid::QCD::TComplex Tp = Grid::sum(Traced);
|
||||
Grid::TComplex Tp = Grid::sum(Traced);
|
||||
Grid::Complex p = Grid::TensorRemove(Tp);
|
||||
|
||||
std::cout<< "Plaq from PlaqStaple "<<p<<std::endl;
|
||||
|
@ -52,8 +52,8 @@ enum ChromaAction {
|
||||
HtContFracZolo
|
||||
};
|
||||
|
||||
void calc_grid (ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag);
|
||||
void calc_chroma (ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag);
|
||||
void calc_grid (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
|
||||
void calc_chroma (ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag);
|
||||
|
||||
#include <chroma.h>
|
||||
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
|
||||
@ -71,10 +71,10 @@ public:
|
||||
typedef LatticeFermion T4;
|
||||
typedef multi1d<LatticeFermion> T5;
|
||||
|
||||
static void ImportGauge(Grid::QCD::LatticeGaugeField & gr,
|
||||
static void ImportGauge(Grid::LatticeGaugeField & gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -112,10 +112,10 @@ public:
|
||||
}}}}
|
||||
}
|
||||
|
||||
static void ImportFermion(Grid::QCD::LatticeFermion & gr,
|
||||
static void ImportFermion(Grid::LatticeFermion & gr,
|
||||
QDP::multi1d<QDP::LatticeFermion> & ch )
|
||||
{
|
||||
Grid::QCD::SpinColourVector F;
|
||||
Grid::SpinColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
QDP::Fermion cF;
|
||||
@ -154,10 +154,10 @@ public:
|
||||
QDP::pokeSite(ch[s],cF,cx);
|
||||
}}}}}
|
||||
}
|
||||
static void ExportFermion(Grid::QCD::LatticeFermion & gr,
|
||||
static void ExportFermion(Grid::LatticeFermion & gr,
|
||||
QDP::multi1d<QDP::LatticeFermion> & ch )
|
||||
{
|
||||
Grid::QCD::SpinColourVector F;
|
||||
Grid::SpinColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
QDP::Fermion cF;
|
||||
@ -384,8 +384,8 @@ int main (int argc,char **argv )
|
||||
* Setup Grid
|
||||
*********************************************************/
|
||||
Grid::Grid_init(&argc,&argv);
|
||||
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridDefaultMpi());
|
||||
|
||||
std::vector<int> gd = UGrid->GlobalDimensions();
|
||||
@ -395,11 +395,11 @@ int main (int argc,char **argv )
|
||||
QDP::Layout::setLattSize(nrow);
|
||||
QDP::Layout::create();
|
||||
|
||||
Grid::GridCartesian * FGrid = Grid::QCD::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
Grid::QCD::LatticeGaugeField lat(UGrid);
|
||||
Grid::QCD::LatticeFermion src(FGrid);
|
||||
Grid::QCD::LatticeFermion res_chroma(FGrid);
|
||||
Grid::QCD::LatticeFermion res_grid (FGrid);
|
||||
Grid::GridCartesian * FGrid = Grid::SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
Grid::LatticeGaugeField lat(UGrid);
|
||||
Grid::LatticeFermion src(FGrid);
|
||||
Grid::LatticeFermion res_chroma(FGrid);
|
||||
Grid::LatticeFermion res_grid (FGrid);
|
||||
|
||||
std::vector<ChromaAction> ActionList({
|
||||
HtCayleyTanh, // Plain old DWF.
|
||||
@ -446,7 +446,7 @@ int main (int argc,char **argv )
|
||||
Chroma::finalize();
|
||||
}
|
||||
|
||||
void calc_chroma(ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag)
|
||||
void calc_chroma(ChromaAction action,Grid::LatticeGaugeField & lat, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag)
|
||||
{
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> u(4);
|
||||
|
||||
@ -483,7 +483,7 @@ void calc_chroma(ChromaAction action,Grid::QCD::LatticeGaugeField & lat, Grid::Q
|
||||
|
||||
|
||||
|
||||
void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res,int dag)
|
||||
void calc_grid(ChromaAction action,Grid::LatticeGaugeField & Umu, Grid::LatticeFermion &src, Grid::LatticeFermion &res,int dag)
|
||||
{
|
||||
using namespace Grid;
|
||||
;
|
||||
@ -493,8 +493,8 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
Grid::GridCartesian * FGrid = (Grid::GridCartesian *) src.Grid();
|
||||
Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
Grid::GridRedBlackCartesian * FrbGrid = Grid::QCD::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
Grid::GridRedBlackCartesian * UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
Grid::GridRedBlackCartesian * FrbGrid = Grid::SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
Grid::GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
@ -502,12 +502,12 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
Grid::gaussian(RNG5,src);
|
||||
Grid::gaussian(RNG5,res);
|
||||
|
||||
Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
|
||||
Grid::SU3::HotConfiguration(RNG4,Umu);
|
||||
|
||||
/*
|
||||
Grid::QCD::LatticeColourMatrix U(UGrid);
|
||||
Grid::LatticeColourMatrix U(UGrid);
|
||||
U=Grid::zero;
|
||||
for(int nn=0;nn<Grid::QCD::Nd;nn++){
|
||||
for(int nn=0;nn<Grid::Nd;nn++){
|
||||
if ( nn>=4 ) {
|
||||
Grid::PokeIndex<LorentzIndex>(Umu,U,nn);
|
||||
}
|
||||
@ -519,7 +519,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
if ( action == HtCayleyTanh ) {
|
||||
|
||||
Grid::QCD::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
|
||||
Grid::DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling domain wall multiply "<<std::endl;
|
||||
|
||||
@ -535,7 +535,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
Grid::Real _b = 0.5*(mobius_scale +1.0);
|
||||
Grid::Real _c = 0.5*(mobius_scale -1.0);
|
||||
Grid::QCD::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
|
||||
Grid::MobiusZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c,zolo_lo,zolo_hi);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling mobius zolo multiply "<<std::endl;
|
||||
|
||||
@ -549,7 +549,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
if ( action == HtCayleyZolo ) {
|
||||
|
||||
Grid::QCD::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
Grid::ShamirZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling shamir zolo multiply "<<std::endl;
|
||||
|
||||
@ -565,7 +565,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
if ( action == HmCayleyTanh ) {
|
||||
Grid::Real _b = 0.5*(mobius_scale +1.0);
|
||||
Grid::Real _c = 0.5*(mobius_scale -1.0);
|
||||
Grid::QCD::MobiusFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c);
|
||||
Grid::MobiusFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,_b,_c);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling mobius tanh multiply "<<std::endl;
|
||||
|
||||
@ -581,7 +581,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
if ( action == HmCayleyTanh ) {
|
||||
|
||||
Grid::QCD::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
|
||||
Grid::ScaledShamirFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,mobius_scale);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling scaled shamir multiply "<<std::endl;
|
||||
|
||||
@ -595,7 +595,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
if ( action == HwCayleyTanh ) {
|
||||
|
||||
Grid::QCD::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
Grid::OverlapWilsonCayleyTanhFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,1.0);
|
||||
|
||||
if ( dag )
|
||||
D.Mdag(src,res);
|
||||
@ -607,7 +607,7 @@ void calc_grid(ChromaAction action,Grid::QCD::LatticeGaugeField & Umu, Grid::QCD
|
||||
|
||||
if ( action == HwCayleyZolo ) {
|
||||
|
||||
Grid::QCD::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
Grid::OverlapWilsonCayleyZolotarevFermionR D(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,_mass,_M5,zolo_lo,zolo_hi);
|
||||
|
||||
if ( dag )
|
||||
D.Mdag(src,res);
|
||||
|
@ -30,8 +30,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
double mq=0.1;
|
||||
|
||||
typedef Grid::QCD::StaggeredImplR::FermionField FermionField;
|
||||
typedef Grid::QCD::LatticeGaugeField GaugeField;
|
||||
typedef Grid::StaggeredImplR::FermionField FermionField;
|
||||
typedef Grid::LatticeGaugeField GaugeField;
|
||||
|
||||
void make_gauge (GaugeField & lat, FermionField &src);
|
||||
void calc_grid (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
|
||||
@ -53,7 +53,7 @@ public:
|
||||
static void ImportGauge(GaugeField & gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -88,7 +88,7 @@ public:
|
||||
static void ExportGauge(GaugeField & gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -124,7 +124,7 @@ public:
|
||||
static void ImportFermion(FermionField & gr,
|
||||
QDP::LatticeStaggeredFermion & ch )
|
||||
{
|
||||
Grid::QCD::ColourVector F;
|
||||
Grid::ColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
|
||||
@ -157,7 +157,7 @@ public:
|
||||
static void ExportFermion(FermionField & gr,
|
||||
QDP::LatticeStaggeredFermion & ch )
|
||||
{
|
||||
Grid::QCD::ColourVector F;
|
||||
Grid::ColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
std::vector<int> x(5);
|
||||
@ -222,8 +222,8 @@ int main (int argc,char **argv )
|
||||
* Setup Grid
|
||||
*********************************************************/
|
||||
Grid::Grid_init(&argc,&argv);
|
||||
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridCartesian * UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::Nd,Grid::vComplex::Nsimd()),
|
||||
Grid::GridDefaultMpi());
|
||||
|
||||
std::vector<int> gd = UGrid->GlobalDimensions();
|
||||
@ -333,7 +333,7 @@ void make_gauge(GaugeField & Umu,FermionField &src)
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu.Grid();
|
||||
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
|
||||
Grid::SU3::HotConfiguration(RNG4,Umu);
|
||||
Grid::gaussian(RNG4,src);
|
||||
}
|
||||
|
||||
@ -343,9 +343,9 @@ void calc_grid(GaugeField & Uthin, GaugeField & Utriple, GaugeField & Ufat, Ferm
|
||||
;
|
||||
|
||||
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Uthin.Grid();
|
||||
Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
Grid::GridRedBlackCartesian * UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
Grid::QCD::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0);
|
||||
Grid::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0);
|
||||
|
||||
std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<<std::endl;
|
||||
|
||||
|
@ -35,8 +35,8 @@
|
||||
double mq = 0.1;
|
||||
|
||||
// Define Wilson Types
|
||||
typedef Grid::QCD::WilsonImplR::FermionField FermionField;
|
||||
typedef Grid::QCD::LatticeGaugeField GaugeField;
|
||||
typedef Grid::WilsonImplR::FermionField FermionField;
|
||||
typedef Grid::LatticeGaugeField GaugeField;
|
||||
|
||||
enum ChromaAction
|
||||
{
|
||||
@ -56,7 +56,7 @@ public:
|
||||
static void ImportGauge(GaugeField &gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -101,7 +101,7 @@ public:
|
||||
static void ExportGauge(GaugeField &gr,
|
||||
QDP::multi1d<QDP::LatticeColorMatrix> &ch)
|
||||
{
|
||||
Grid::QCD::LorentzColourMatrix LCM;
|
||||
Grid::LorentzColourMatrix LCM;
|
||||
Grid::Complex cc;
|
||||
QDP::ColorMatrix cm;
|
||||
QDP::Complex c;
|
||||
@ -145,10 +145,10 @@ public:
|
||||
}
|
||||
|
||||
// Specific for Wilson Fermions
|
||||
static void ImportFermion(Grid::QCD::LatticeFermion &gr,
|
||||
static void ImportFermion(Grid::LatticeFermion &gr,
|
||||
QDP::LatticeFermion &ch)
|
||||
{
|
||||
Grid::QCD::SpinColourVector F;
|
||||
Grid::SpinColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
QDP::Fermion cF;
|
||||
@ -195,10 +195,10 @@ public:
|
||||
}
|
||||
|
||||
// Specific for 4d Wilson fermions
|
||||
static void ExportFermion(Grid::QCD::LatticeFermion &gr,
|
||||
static void ExportFermion(Grid::LatticeFermion &gr,
|
||||
QDP::LatticeFermion &ch)
|
||||
{
|
||||
Grid::QCD::SpinColourVector F;
|
||||
Grid::SpinColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
QDP::Fermion cF;
|
||||
@ -342,19 +342,18 @@ void calc_chroma(ChromaAction action, GaugeField &lat, FermionField &src, Fermio
|
||||
void make_gauge(GaugeField &Umu, FermionField &src)
|
||||
{
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
std::vector<int> seeds4({1, 2, 3, 4});
|
||||
|
||||
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
|
||||
Grid::GridParallelRNG RNG4(UGrid);
|
||||
RNG4.SeedFixedIntegers(seeds4);
|
||||
Grid::QCD::SU3::HotConfiguration(RNG4, Umu);
|
||||
Grid::SU3::HotConfiguration(RNG4, Umu);
|
||||
|
||||
// Fermion field
|
||||
Grid::gaussian(RNG4, src);
|
||||
/*
|
||||
Grid::QCD::SpinColourVector F;
|
||||
Grid::SpinColourVector F;
|
||||
Grid::Complex c;
|
||||
|
||||
|
||||
@ -391,13 +390,12 @@ void make_gauge(GaugeField &Umu, FermionField &src)
|
||||
*/
|
||||
}
|
||||
|
||||
void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD::LatticeFermion &src, Grid::QCD::LatticeFermion &res, int dag)
|
||||
void calc_grid(ChromaAction action, Grid::LatticeGaugeField &Umu, Grid::LatticeFermion &src, Grid::LatticeFermion &res, int dag)
|
||||
{
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
Grid::GridCartesian *UGrid = (Grid::GridCartesian *)Umu._grid;
|
||||
Grid::GridRedBlackCartesian *UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
Grid::GridRedBlackCartesian *UrbGrid = Grid::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
Grid::RealD _mass = mq;
|
||||
|
||||
@ -409,7 +407,7 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD
|
||||
anis.xi_0 = 2.0;
|
||||
anis.nu = 1.0;
|
||||
WilsonImplParams iParam;
|
||||
Grid::QCD::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis);
|
||||
Grid::WilsonFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, iParam, anis);
|
||||
|
||||
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Fermion multiply " << std::endl;
|
||||
|
||||
@ -430,7 +428,7 @@ void calc_grid(ChromaAction action, Grid::QCD::LatticeGaugeField &Umu, Grid::QCD
|
||||
anis.xi_0 = 2.0;
|
||||
anis.nu = 1.0;
|
||||
WilsonImplParams CloverImplParam;
|
||||
Grid::QCD::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam);
|
||||
Grid::WilsonCloverFermionR Wf(Umu, *UGrid, *UrbGrid, _mass, _csw_r, _csw_t, anis, CloverImplParam);
|
||||
Wf.ImportGauge(Umu);
|
||||
|
||||
std::cout << Grid::GridLogMessage << " Calling Grid Wilson Clover Fermion multiply " << std::endl;
|
||||
@ -458,8 +456,8 @@ int main(int argc, char **argv)
|
||||
* Setup Grid
|
||||
*********************************************************/
|
||||
Grid::Grid_init(&argc, &argv);
|
||||
Grid::GridCartesian *UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::QCD::Nd, Grid::vComplex::Nsimd()),
|
||||
Grid::GridCartesian *UGrid = Grid::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
|
||||
Grid::GridDefaultSimd(Grid::Nd, Grid::vComplex::Nsimd()),
|
||||
Grid::GridDefaultMpi());
|
||||
|
||||
std::vector<int> gd = UGrid->GlobalDimensions();
|
||||
|
@ -30,8 +30,6 @@ Author: Peter Boyle <paboyle@ph.ed.ac.uk>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
|
@ -32,7 +32,6 @@
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
int main(int argc, char **argv) {
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user