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

Faster derivative for WilsonGauge

Browse Source
This commit is contained in:
Guido Cossu 2017-05-26 14:31:49 +01:00
parent ffb91e53d2
commit 0de314870d
2 changed files with 37 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
lib/qcd/action/gauge/WilsonGaugeAction.h
View File

@ -71,14 +71,18 @@ class WilsonGaugeAction : public Action<typename Gimpl::GaugeField> {
RealD factor = 0.5 * beta / RealD(Nc); RealD factor = 0.5 * beta / RealD(Nc);
GaugeLinkField Umu(U._grid); //GaugeLinkField Umu(U._grid);
GaugeLinkField dSdU_mu(U._grid); GaugeLinkField dSdU_mu(U._grid);
for (int mu = 0; mu < Nd; mu++) { for (int mu = 0; mu < Nd; mu++) {
Umu = PeekIndex<LorentzIndex>(U, mu); //Umu = PeekIndex<LorentzIndex>(U, mu);
// Staple in direction mu // Staple in direction mu
WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu); //WilsonLoops<Gimpl>::Staple(dSdU_mu, U, mu);
dSdU_mu = Ta(Umu * dSdU_mu) * factor; //dSdU_mu = Ta(Umu * dSdU_mu) * factor;
WilsonLoops<Gimpl>::StapleMult(dSdU_mu, U, mu);
dSdU_mu = Ta(dSdU_mu) * factor;
PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu); PokeIndex<LorentzIndex>(dSdU, dSdU_mu, mu);
} }

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

@ -188,13 +188,10 @@ public:
} }
} }
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
// For the force term
static void StapleMult(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
GridBase *grid = Umu._grid; GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd, grid); std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) { for (int d = 0; d < Nd; d++) {
// this operation is taking too much time // this operation is taking too much time
@ -204,6 +201,31 @@ public:
GaugeMat tmp1(grid); GaugeMat tmp1(grid);
GaugeMat tmp2(grid); GaugeMat tmp2(grid);
for (int nu = 0; nu < Nd; nu++) {
if (nu != mu) {
// this is ~10% faster than the Staple
tmp1 = Cshift(U[nu], mu, 1);
tmp2 = Cshift(U[mu], nu, 1);
staple += tmp1* adj(U[nu]*tmp2);
tmp2 = adj(U[mu]*tmp1)*U[nu];
staple += Cshift(tmp2, nu, -1);
}
}
staple = U[mu]*staple;
}
//////////////////////////////////////////////////
// the sum over all staples on each site
//////////////////////////////////////////////////
static void Staple(GaugeMat &staple, const GaugeLorentz &Umu, int mu) {
GridBase *grid = Umu._grid;
std::vector<GaugeMat> U(Nd, grid);
for (int d = 0; d < Nd; d++) {
U[d] = PeekIndex<LorentzIndex>(Umu, d);
}
staple = zero;
for (int nu = 0; nu < Nd; nu++) { for (int nu = 0; nu < Nd; nu++) {
@ -217,34 +239,23 @@ public:
// | // |
// __| // __|
// //
tmp1 = Cshift(U[nu], mu, 1);
tmp2 = Cshift(U[mu], nu, 1);
staple += tmp1*adj(U[nu] * tmp2);
/*
staple += Gimpl::ShiftStaple( staple += Gimpl::ShiftStaple(
Gimpl::CovShiftForward( Gimpl::CovShiftForward(
U[nu], nu, U[nu], nu,
Gimpl::CovShiftBackward( Gimpl::CovShiftBackward(
U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))), U[mu], mu, Gimpl::CovShiftIdentityBackward(U[nu], nu))),
mu); mu);
*/
// __ // __
// | // |
// |__ // |__
// //
// //
tmp2 = adj(U[mu]*tmp1)*U[nu];
staple += Cshift(tmp2, nu, -1);
/*
staple += Gimpl::ShiftStaple( staple += Gimpl::ShiftStaple(
Gimpl::CovShiftBackward(U[nu], nu, Gimpl::CovShiftBackward(U[nu], nu,
Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu); Gimpl::CovShiftBackward(U[mu], mu, U[nu])), mu);
*/
} }
} }
} }