From f77f3a6598decb2c4e3ef55bcedd1dacbc3623d0 Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Wed, 6 Apr 2022 10:21:04 -0400 Subject: [PATCH 1/4] Imported G-parity flavor algebra + tester from feature/gparity_HMC branch --- Grid/GridQCDcore.h | 1 + Grid/qcd/QCD.h | 25 ++ Grid/qcd/gparity/Gparity.h | 6 + Grid/qcd/gparity/GparityFlavour.cc | 34 +++ Grid/qcd/gparity/GparityFlavour.h | 475 +++++++++++++++++++++++++++++ tests/core/Test_gparity_flavour.cc | 177 +++++++++++ 6 files changed, 718 insertions(+) create mode 100644 Grid/qcd/gparity/Gparity.h create mode 100644 Grid/qcd/gparity/GparityFlavour.cc create mode 100644 Grid/qcd/gparity/GparityFlavour.h create mode 100644 tests/core/Test_gparity_flavour.cc diff --git a/Grid/GridQCDcore.h b/Grid/GridQCDcore.h index cae6f43f..065b62cd 100644 --- a/Grid/GridQCDcore.h +++ b/Grid/GridQCDcore.h @@ -36,6 +36,7 @@ Author: paboyle #include #include #include +#include #include #include NAMESPACE_CHECK(GridQCDCore); diff --git a/Grid/qcd/QCD.h b/Grid/qcd/QCD.h index 858aead7..81356a66 100644 --- a/Grid/qcd/QCD.h +++ b/Grid/qcd/QCD.h @@ -63,6 +63,7 @@ static constexpr int Ngp=2; // gparity index range #define ColourIndex (2) #define SpinIndex (1) #define LorentzIndex (0) +#define GparityFlavourIndex (0) // Also should make these a named enum type static constexpr int DaggerNo=0; @@ -87,6 +88,8 @@ template struct isCoarsened { template using IfCoarsened = Invoke::value,int> > ; template using IfNotCoarsened = Invoke::value,int> > ; +const int GparityFlavourTensorIndex = 3; //TensorLevel counts from the bottom! + // ChrisK very keen to add extra space for Gparity doubling. // // Also add domain wall index, in a way where Wilson operator @@ -110,8 +113,10 @@ template using iHalfSpinColourVector = iScalar using iSpinColourSpinColourMatrix = iScalar, Ns>, Nc>, Ns> >; +template using iGparityFlavourVector = iVector >, Ngp>; template using iGparitySpinColourVector = iVector, Ns>, Ngp >; template using iGparityHalfSpinColourVector = iVector, Nhs>, Ngp >; +template using iGparityFlavourMatrix = iMatrix >, Ngp>; // Spin matrix typedef iSpinMatrix SpinMatrix; @@ -176,6 +181,16 @@ typedef iDoubleStoredColourMatrix vDoubleStoredColourMatrix; typedef iDoubleStoredColourMatrix vDoubleStoredColourMatrixF; typedef iDoubleStoredColourMatrix vDoubleStoredColourMatrixD; +//G-parity flavour matrix +typedef iGparityFlavourMatrix GparityFlavourMatrix; +typedef iGparityFlavourMatrix GparityFlavourMatrixF; +typedef iGparityFlavourMatrix GparityFlavourMatrixD; + +typedef iGparityFlavourMatrix vGparityFlavourMatrix; +typedef iGparityFlavourMatrix vGparityFlavourMatrixF; +typedef iGparityFlavourMatrix vGparityFlavourMatrixD; + + // Spin vector typedef iSpinVector SpinVector; typedef iSpinVector SpinVectorF; @@ -220,6 +235,16 @@ typedef iHalfSpinColourVector HalfSpinColourVectorD; typedef iHalfSpinColourVector vHalfSpinColourVector; typedef iHalfSpinColourVector vHalfSpinColourVectorF; typedef iHalfSpinColourVector vHalfSpinColourVectorD; + +//G-parity flavour vector +typedef iGparityFlavourVector GparityFlavourVector; +typedef iGparityFlavourVector GparityFlavourVectorF; +typedef iGparityFlavourVector GparityFlavourVectorD; + +typedef iGparityFlavourVector vGparityFlavourVector; +typedef iGparityFlavourVector vGparityFlavourVectorF; +typedef iGparityFlavourVector vGparityFlavourVectorD; + // singlets typedef iSinglet TComplex; // FIXME This is painful. Tensor singlet complex type. diff --git a/Grid/qcd/gparity/Gparity.h b/Grid/qcd/gparity/Gparity.h new file mode 100644 index 00000000..ce1c70eb --- /dev/null +++ b/Grid/qcd/gparity/Gparity.h @@ -0,0 +1,6 @@ +#ifndef GRID_GPARITY_H_ +#define GRID_GPARITY_H_ + +#include + +#endif diff --git a/Grid/qcd/gparity/GparityFlavour.cc b/Grid/qcd/gparity/GparityFlavour.cc new file mode 100644 index 00000000..4596f96b --- /dev/null +++ b/Grid/qcd/gparity/GparityFlavour.cc @@ -0,0 +1,34 @@ +#include + +NAMESPACE_BEGIN(Grid); + +const std::array GparityFlavour::sigma_mu = {{ + GparityFlavour(GparityFlavour::Algebra::SigmaX), + GparityFlavour(GparityFlavour::Algebra::SigmaY), + GparityFlavour(GparityFlavour::Algebra::SigmaZ) + }}; + +const std::array GparityFlavour::sigma_all = {{ + GparityFlavour(GparityFlavour::Algebra::Identity), + GparityFlavour(GparityFlavour::Algebra::SigmaX), + GparityFlavour(GparityFlavour::Algebra::SigmaY), + GparityFlavour(GparityFlavour::Algebra::SigmaZ), + GparityFlavour(GparityFlavour::Algebra::ProjPlus), + GparityFlavour(GparityFlavour::Algebra::ProjMinus) +}}; + +const std::array GparityFlavour::name = {{ + "SigmaX", + "MinusSigmaX", + "SigmaY", + "MinusSigmaY", + "SigmaZ", + "MinusSigmaZ", + "Identity", + "MinusIdentity", + "ProjPlus", + "MinusProjPlus", + "ProjMinus", + "MinusProjMinus"}}; + +NAMESPACE_END(Grid); diff --git a/Grid/qcd/gparity/GparityFlavour.h b/Grid/qcd/gparity/GparityFlavour.h new file mode 100644 index 00000000..b2009235 --- /dev/null +++ b/Grid/qcd/gparity/GparityFlavour.h @@ -0,0 +1,475 @@ +#ifndef GRID_QCD_GPARITY_FLAVOUR_H +#define GRID_QCD_GPARITY_FLAVOUR_H + +//Support for flavour-matrix operations acting on the G-parity flavour index + +#include + +NAMESPACE_BEGIN(Grid); + +class GparityFlavour { + public: + GRID_SERIALIZABLE_ENUM(Algebra, undef, + SigmaX, 0, + MinusSigmaX, 1, + SigmaY, 2, + MinusSigmaY, 3, + SigmaZ, 4, + MinusSigmaZ, 5, + Identity, 6, + MinusIdentity, 7, + ProjPlus, 8, + MinusProjPlus, 9, + ProjMinus, 10, + MinusProjMinus, 11 + ); + static constexpr unsigned int nSigma = 12; + static const std::array name; + static const std::array sigma_mu; + static const std::array sigma_all; + Algebra g; + public: + accelerator GparityFlavour(Algebra initg): g(initg) {} +}; + + + +// 0 1 x vector +// 1 0 +template +accelerator_inline void multFlavourSigmaX(iVector &ret, const iVector &rhs) +{ + ret(0) = rhs(1); + ret(1) = rhs(0); +}; +template +accelerator_inline void lmultFlavourSigmaX(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(1,0); + ret(0,1) = rhs(1,1); + ret(1,0) = rhs(0,0); + ret(1,1) = rhs(0,1); +}; +template +accelerator_inline void rmultFlavourSigmaX(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(0,1); + ret(0,1) = rhs(0,0); + ret(1,0) = rhs(1,1); + ret(1,1) = rhs(1,0); +}; + + +template +accelerator_inline void multFlavourMinusSigmaX(iVector &ret, const iVector &rhs) +{ + ret(0) = -rhs(1); + ret(1) = -rhs(0); +}; +template +accelerator_inline void lmultFlavourMinusSigmaX(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(1,0); + ret(0,1) = -rhs(1,1); + ret(1,0) = -rhs(0,0); + ret(1,1) = -rhs(0,1); +}; +template +accelerator_inline void rmultFlavourMinusSigmaX(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(0,1); + ret(0,1) = -rhs(0,0); + ret(1,0) = -rhs(1,1); + ret(1,1) = -rhs(1,0); +}; + + + + + +// 0 -i x vector +// i 0 +template +accelerator_inline void multFlavourSigmaY(iVector &ret, const iVector &rhs) +{ + ret(0) = timesMinusI(rhs(1)); + ret(1) = timesI(rhs(0)); +}; +template +accelerator_inline void lmultFlavourSigmaY(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = timesMinusI(rhs(1,0)); + ret(0,1) = timesMinusI(rhs(1,1)); + ret(1,0) = timesI(rhs(0,0)); + ret(1,1) = timesI(rhs(0,1)); +}; +template +accelerator_inline void rmultFlavourSigmaY(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = timesI(rhs(0,1)); + ret(0,1) = timesMinusI(rhs(0,0)); + ret(1,0) = timesI(rhs(1,1)); + ret(1,1) = timesMinusI(rhs(1,0)); +}; + +template +accelerator_inline void multFlavourMinusSigmaY(iVector &ret, const iVector &rhs) +{ + ret(0) = timesI(rhs(1)); + ret(1) = timesMinusI(rhs(0)); +}; +template +accelerator_inline void lmultFlavourMinusSigmaY(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = timesI(rhs(1,0)); + ret(0,1) = timesI(rhs(1,1)); + ret(1,0) = timesMinusI(rhs(0,0)); + ret(1,1) = timesMinusI(rhs(0,1)); +}; +template +accelerator_inline void rmultFlavourMinusSigmaY(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = timesMinusI(rhs(0,1)); + ret(0,1) = timesI(rhs(0,0)); + ret(1,0) = timesMinusI(rhs(1,1)); + ret(1,1) = timesI(rhs(1,0)); +}; + + + + + +// 1 0 x vector +// 0 -1 +template +accelerator_inline void multFlavourSigmaZ(iVector &ret, const iVector &rhs) +{ + ret(0) = rhs(0); + ret(1) = -rhs(1); +}; +template +accelerator_inline void lmultFlavourSigmaZ(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(0,0); + ret(0,1) = rhs(0,1); + ret(1,0) = -rhs(1,0); + ret(1,1) = -rhs(1,1); +}; +template +accelerator_inline void rmultFlavourSigmaZ(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(0,0); + ret(0,1) = -rhs(0,1); + ret(1,0) = rhs(1,0); + ret(1,1) = -rhs(1,1); +}; + + +template +accelerator_inline void multFlavourMinusSigmaZ(iVector &ret, const iVector &rhs) +{ + ret(0) = -rhs(0); + ret(1) = rhs(1); +}; +template +accelerator_inline void lmultFlavourMinusSigmaZ(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(0,0); + ret(0,1) = -rhs(0,1); + ret(1,0) = rhs(1,0); + ret(1,1) = rhs(1,1); +}; +template +accelerator_inline void rmultFlavourMinusSigmaZ(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(0,0); + ret(0,1) = rhs(0,1); + ret(1,0) = -rhs(1,0); + ret(1,1) = rhs(1,1); +}; + + + + + + +template +accelerator_inline void multFlavourIdentity(iVector &ret, const iVector &rhs) +{ + ret(0) = rhs(0); + ret(1) = rhs(1); +}; +template +accelerator_inline void lmultFlavourIdentity(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(0,0); + ret(0,1) = rhs(0,1); + ret(1,0) = rhs(1,0); + ret(1,1) = rhs(1,1); +}; +template +accelerator_inline void rmultFlavourIdentity(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = rhs(0,0); + ret(0,1) = rhs(0,1); + ret(1,0) = rhs(1,0); + ret(1,1) = rhs(1,1); +}; + +template +accelerator_inline void multFlavourMinusIdentity(iVector &ret, const iVector &rhs) +{ + ret(0) = -rhs(0); + ret(1) = -rhs(1); +}; +template +accelerator_inline void lmultFlavourMinusIdentity(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(0,0); + ret(0,1) = -rhs(0,1); + ret(1,0) = -rhs(1,0); + ret(1,1) = -rhs(1,1); +}; +template +accelerator_inline void rmultFlavourMinusIdentity(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -rhs(0,0); + ret(0,1) = -rhs(0,1); + ret(1,0) = -rhs(1,0); + ret(1,1) = -rhs(1,1); +}; + + + + + +//G-parity flavour projection 1/2(1+\sigma_2) +//1 -i +//i 1 +template +accelerator_inline void multFlavourProjPlus(iVector &ret, const iVector &rhs) +{ + ret(0) = 0.5*rhs(0) + 0.5*timesMinusI(rhs(1)); + ret(1) = 0.5*timesI(rhs(0)) + 0.5*rhs(1); +}; +template +accelerator_inline void lmultFlavourProjPlus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = 0.5*rhs(0,0) + 0.5*timesMinusI(rhs(1,0)); + ret(0,1) = 0.5*rhs(0,1) + 0.5*timesMinusI(rhs(1,1)); + ret(1,0) = 0.5*timesI(rhs(0,0)) + 0.5*rhs(1,0); + ret(1,1) = 0.5*timesI(rhs(0,1)) + 0.5*rhs(1,1); +}; +template +accelerator_inline void rmultFlavourProjPlus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = 0.5*rhs(0,0) + 0.5*timesI(rhs(0,1)); + ret(0,1) = 0.5*timesMinusI(rhs(0,0)) + 0.5*rhs(0,1); + ret(1,0) = 0.5*rhs(1,0) + 0.5*timesI(rhs(1,1)); + ret(1,1) = 0.5*timesMinusI(rhs(1,0)) + 0.5*rhs(1,1); +}; + + +template +accelerator_inline void multFlavourMinusProjPlus(iVector &ret, const iVector &rhs) +{ + ret(0) = -0.5*rhs(0) + 0.5*timesI(rhs(1)); + ret(1) = 0.5*timesMinusI(rhs(0)) - 0.5*rhs(1); +}; +template +accelerator_inline void lmultFlavourMinusProjPlus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -0.5*rhs(0,0) + 0.5*timesI(rhs(1,0)); + ret(0,1) = -0.5*rhs(0,1) + 0.5*timesI(rhs(1,1)); + ret(1,0) = 0.5*timesMinusI(rhs(0,0)) - 0.5*rhs(1,0); + ret(1,1) = 0.5*timesMinusI(rhs(0,1)) - 0.5*rhs(1,1); +}; +template +accelerator_inline void rmultFlavourMinusProjPlus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -0.5*rhs(0,0) + 0.5*timesMinusI(rhs(0,1)); + ret(0,1) = 0.5*timesI(rhs(0,0)) - 0.5*rhs(0,1); + ret(1,0) = -0.5*rhs(1,0) + 0.5*timesMinusI(rhs(1,1)); + ret(1,1) = 0.5*timesI(rhs(1,0)) - 0.5*rhs(1,1); +}; + + + + + +//G-parity flavour projection 1/2(1-\sigma_2) +//1 i +//-i 1 +template +accelerator_inline void multFlavourProjMinus(iVector &ret, const iVector &rhs) +{ + ret(0) = 0.5*rhs(0) + 0.5*timesI(rhs(1)); + ret(1) = 0.5*timesMinusI(rhs(0)) + 0.5*rhs(1); +}; +template +accelerator_inline void lmultFlavourProjMinus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = 0.5*rhs(0,0) + 0.5*timesI(rhs(1,0)); + ret(0,1) = 0.5*rhs(0,1) + 0.5*timesI(rhs(1,1)); + ret(1,0) = 0.5*timesMinusI(rhs(0,0)) + 0.5*rhs(1,0); + ret(1,1) = 0.5*timesMinusI(rhs(0,1)) + 0.5*rhs(1,1); +}; +template +accelerator_inline void rmultFlavourProjMinus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = 0.5*rhs(0,0) + 0.5*timesMinusI(rhs(0,1)); + ret(0,1) = 0.5*timesI(rhs(0,0)) + 0.5*rhs(0,1); + ret(1,0) = 0.5*rhs(1,0) + 0.5*timesMinusI(rhs(1,1)); + ret(1,1) = 0.5*timesI(rhs(1,0)) + 0.5*rhs(1,1); +}; + + +template +accelerator_inline void multFlavourMinusProjMinus(iVector &ret, const iVector &rhs) +{ + ret(0) = -0.5*rhs(0) + 0.5*timesMinusI(rhs(1)); + ret(1) = 0.5*timesI(rhs(0)) - 0.5*rhs(1); +}; +template +accelerator_inline void lmultFlavourMinusProjMinus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -0.5*rhs(0,0) + 0.5*timesMinusI(rhs(1,0)); + ret(0,1) = -0.5*rhs(0,1) + 0.5*timesMinusI(rhs(1,1)); + ret(1,0) = 0.5*timesI(rhs(0,0)) - 0.5*rhs(1,0); + ret(1,1) = 0.5*timesI(rhs(0,1)) - 0.5*rhs(1,1); +}; +template +accelerator_inline void rmultFlavourMinusProjMinus(iMatrix &ret, const iMatrix &rhs) +{ + ret(0,0) = -0.5*rhs(0,0) + 0.5*timesI(rhs(0,1)); + ret(0,1) = 0.5*timesMinusI(rhs(0,0)) - 0.5*rhs(0,1); + ret(1,0) = -0.5*rhs(1,0) + 0.5*timesI(rhs(1,1)); + ret(1,1) = 0.5*timesMinusI(rhs(1,0)) - 0.5*rhs(1,1); +}; + + + + + + + + + + +template +accelerator_inline auto operator*(const GparityFlavour &G, const iVector &arg) +->typename std::enable_if, GparityFlavourTensorIndex>::value, iVector>::type +{ + iVector ret; + + switch (G.g) + { + case GparityFlavour::Algebra::SigmaX: + multFlavourSigmaX(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaX: + multFlavourMinusSigmaX(ret, arg); break; + case GparityFlavour::Algebra::SigmaY: + multFlavourSigmaY(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaY: + multFlavourMinusSigmaY(ret, arg); break; + case GparityFlavour::Algebra::SigmaZ: + multFlavourSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaZ: + multFlavourMinusSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::Identity: + multFlavourIdentity(ret, arg); break; + case GparityFlavour::Algebra::MinusIdentity: + multFlavourMinusIdentity(ret, arg); break; + case GparityFlavour::Algebra::ProjPlus: + multFlavourProjPlus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjPlus: + multFlavourMinusProjPlus(ret, arg); break; + case GparityFlavour::Algebra::ProjMinus: + multFlavourProjMinus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjMinus: + multFlavourMinusProjMinus(ret, arg); break; + default: assert(0); + } + + return ret; +} + +template +accelerator_inline auto operator*(const GparityFlavour &G, const iMatrix &arg) +->typename std::enable_if, GparityFlavourTensorIndex>::value, iMatrix>::type +{ + iMatrix ret; + + switch (G.g) + { + case GparityFlavour::Algebra::SigmaX: + lmultFlavourSigmaX(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaX: + lmultFlavourMinusSigmaX(ret, arg); break; + case GparityFlavour::Algebra::SigmaY: + lmultFlavourSigmaY(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaY: + lmultFlavourMinusSigmaY(ret, arg); break; + case GparityFlavour::Algebra::SigmaZ: + lmultFlavourSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaZ: + lmultFlavourMinusSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::Identity: + lmultFlavourIdentity(ret, arg); break; + case GparityFlavour::Algebra::MinusIdentity: + lmultFlavourMinusIdentity(ret, arg); break; + case GparityFlavour::Algebra::ProjPlus: + lmultFlavourProjPlus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjPlus: + lmultFlavourMinusProjPlus(ret, arg); break; + case GparityFlavour::Algebra::ProjMinus: + lmultFlavourProjMinus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjMinus: + lmultFlavourMinusProjMinus(ret, arg); break; + default: assert(0); + } + + return ret; +} + +template +accelerator_inline auto operator*(const iMatrix &arg, const GparityFlavour &G) +->typename std::enable_if, GparityFlavourTensorIndex>::value, iMatrix>::type +{ + iMatrix ret; + + switch (G.g) + { + case GparityFlavour::Algebra::SigmaX: + rmultFlavourSigmaX(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaX: + rmultFlavourMinusSigmaX(ret, arg); break; + case GparityFlavour::Algebra::SigmaY: + rmultFlavourSigmaY(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaY: + rmultFlavourMinusSigmaY(ret, arg); break; + case GparityFlavour::Algebra::SigmaZ: + rmultFlavourSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::MinusSigmaZ: + rmultFlavourMinusSigmaZ(ret, arg); break; + case GparityFlavour::Algebra::Identity: + rmultFlavourIdentity(ret, arg); break; + case GparityFlavour::Algebra::MinusIdentity: + rmultFlavourMinusIdentity(ret, arg); break; + case GparityFlavour::Algebra::ProjPlus: + rmultFlavourProjPlus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjPlus: + rmultFlavourMinusProjPlus(ret, arg); break; + case GparityFlavour::Algebra::ProjMinus: + rmultFlavourProjMinus(ret, arg); break; + case GparityFlavour::Algebra::MinusProjMinus: + rmultFlavourMinusProjMinus(ret, arg); break; + default: assert(0); + } + + return ret; +} + +NAMESPACE_END(Grid); + +#endif // include guard diff --git a/tests/core/Test_gparity_flavour.cc b/tests/core/Test_gparity_flavour.cc new file mode 100644 index 00000000..5b204e04 --- /dev/null +++ b/tests/core/Test_gparity_flavour.cc @@ -0,0 +1,177 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: ./tests/Test_gparity_flavour.cc + +Copyright (C) 2015-2017 + +Author: Christopher Kelly +Author: Peter Boyle + +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 +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ + /* END LEGAL */ +#include + +using namespace Grid; + +static constexpr double tolerance = 1.0e-6; +static std::array testAlgebra; + +void print(const GparityFlavourMatrix &g) +{ + for(int i = 0; i < Ngp; i++) + { + std::cout << GridLogMessage << "("; + for(int j=0;j testg; + const Complex I(0., 1.), mI(0., -1.); + + // 0 1 + // 1 0 + testg[0] = Zero(); + testg[0](0, 1)()() = 1.; + testg[0](1, 0)()() = 1.; + std::cout << GridLogMessage << "test SigmaX= " << std::endl; + print(testg[0]); + + // 0 -i + // i 0 + testg[1] = Zero(); + testg[1](0, 1)()() = mI; + testg[1](1, 0)()() = I; + std::cout << GridLogMessage << "test SigmaY= " << std::endl; + print(testg[1]); + + // 1 0 + // 0 -1 + testg[2] = Zero(); + testg[2](0, 0)()() = 1.0; + testg[2](1, 1)()() = -1.0; + std::cout << GridLogMessage << "test SigmaZ= " << std::endl; + print(testg[2]); + + +#define DEFINE_TEST_G(g, exp)\ +testAlgebra[GparityFlavour::Algebra::g] = exp; \ +testAlgebra[GparityFlavour::Algebra::Minus##g] = -exp; + + DEFINE_TEST_G(SigmaX , testg[0]); + DEFINE_TEST_G(SigmaY , testg[1]); + DEFINE_TEST_G(SigmaZ , testg[2]); + DEFINE_TEST_G(Identity , 1.); + + GparityFlavourMatrix pplus; + pplus = 1.0; + pplus = pplus + testg[1]; + pplus = pplus * 0.5; + + DEFINE_TEST_G(ProjPlus , pplus); + + GparityFlavourMatrix pminus; + pminus = 1.0; + pminus = pminus - testg[1]; + pminus = pminus * 0.5; + + DEFINE_TEST_G(ProjMinus , pminus); + +#undef DEFINE_TEST_G +} + +template +void test(const Expr &a, const Expr &b) +{ + if (norm2(a - b) < tolerance) + { + std::cout << "[OK] "; + } + else + { + std::cout << "[fail]" << std::endl; + std::cout << GridLogError << "a= " << a << std::endl; + std::cout << GridLogError << "is different (tolerance= " << tolerance << ") from " << std::endl; + std::cout << GridLogError << "b= " << b << std::endl; + exit(EXIT_FAILURE); + } +} + +void checkSigma(const GparityFlavour::Algebra a, GridSerialRNG &rng) +{ + GparityFlavourVector v; + GparityFlavourMatrix m, &testg = testAlgebra[a]; + GparityFlavour g(a); + + random(rng, v); + random(rng, m); + + std::cout << GridLogMessage << "Checking " << GparityFlavour::name[a] << ": "; + std::cout << "vecmul "; + test(g*v, testg*v); + std::cout << "matlmul "; + test(g*m, testg*m); + std::cout << "matrmul "; + test(m*g, m*testg); + std::cout << std::endl; +} + +int main(int argc, char *argv[]) +{ + Grid_init(&argc,&argv); + + Coordinate latt_size = GridDefaultLatt(); + Coordinate simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); + Coordinate mpi_layout = GridDefaultMpi(); + + GridCartesian Grid(latt_size,simd_layout,mpi_layout); + GridSerialRNG sRNG; + + sRNG.SeedFixedIntegers(std::vector({45,12,81,9})); + + std::cout << GridLogMessage << "======== Test algebra" << std::endl; + createTestAlgebra(); + std::cout << GridLogMessage << "======== Multiplication operators check" << std::endl; + for (int i = 0; i < GparityFlavour::nSigma; ++i) + { + checkSigma(i, sRNG); + } + std::cout << GridLogMessage << std::endl; + + Grid_finalize(); + + return EXIT_SUCCESS; +} From 81fe4c937e61a26c35db1a92ba055b18e15045eb Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Tue, 12 Apr 2022 09:51:59 -0400 Subject: [PATCH 2/4] Hopefully fix link errors on Intel compilers due to having no function body for MomentumFilterBase::apply_phase --- Grid/qcd/action/filters/MomentumFilter.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Grid/qcd/action/filters/MomentumFilter.h b/Grid/qcd/action/filters/MomentumFilter.h index 2a15d80c..864166f5 100644 --- a/Grid/qcd/action/filters/MomentumFilter.h +++ b/Grid/qcd/action/filters/MomentumFilter.h @@ -37,7 +37,7 @@ NAMESPACE_BEGIN(Grid); template struct MomentumFilterBase{ - virtual void applyFilter(MomentaField &P) const; + virtual void applyFilter(MomentaField &P) const = 0; }; //Do nothing From 6121397587adfbaf876ce80dc697b631c022929e Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Mon, 9 May 2022 16:27:57 -0400 Subject: [PATCH 3/4] Imported changes from feature/gparity_HMC branch: Added storage of final true residual in mixed-prec CG and enhanced log output Fixed const correctness of multi-shift constructor Added a mixed precision variant of the multi-shift algorithm that uses a single precision operator and applies periodic reliable update to the residual Added tests/solver/Test_dwf_multishift_mixedprec to test the above Fixed local coherence lanczos using the (large!) max approx to the chebyshev eval as the scale from which to judge the quality of convergence, resulting a test that always passes Added a method to local coherence lanczos class that returns the fine eval/evec pair Added iterative log output to power method Added optional disabling of the plaquette check in Nerscio to support loading old G-parity configs which have a factor of 2 error in the plaquette G-parity Dirac op no longer allows GPBC in the time direction; instead we toggle between periodic and antiperiodic Replaced thread_for G-parity 5D force insertion implementation with accelerator_for version capable of running on GPUs Generalized tests/lanczos/Test_dwf_lanczos to support regular DWF as well as Gparity, with the action chosen by a command line option Modified tests/forces/Test_dwf_gpforce,Test_gpdwf_force,Test_gpwilson_force to use GPBC a spatial direction rather than the t-direction, and antiperiodic BCs for time direction tests/core/Test_gparity now supports using APBC in time direction using command line toggle --- Grid/algorithms/Algorithms.h | 1 + .../iterative/ConjugateGradientMixedPrec.h | 5 + .../iterative/ConjugateGradientMultiShift.h | 5 +- .../ConjugateGradientMultiShiftMixedPrec.h | 409 ++++++++++++++++++ .../iterative/LocalCoherenceLanczos.h | 48 +- Grid/algorithms/iterative/PowerMethod.h | 2 + Grid/parallelIO/NerscIO.h | 6 +- Grid/qcd/action/fermion/GparityWilsonImpl.h | 136 ++++-- tests/core/Test_gparity.cc | 36 +- tests/forces/Test_dwf_gpforce.cc | 18 +- tests/forces/Test_gpdwf_force.cc | 6 +- tests/forces/Test_gpwilson_force.cc | 8 +- tests/lanczos/Test_dwf_lanczos.cc | 81 ++-- tests/solver/Test_dwf_multishift_mixedprec.cc | 184 ++++++++ 14 files changed, 852 insertions(+), 93 deletions(-) create mode 100644 Grid/algorithms/iterative/ConjugateGradientMultiShiftMixedPrec.h create mode 100644 tests/solver/Test_dwf_multishift_mixedprec.cc diff --git a/Grid/algorithms/Algorithms.h b/Grid/algorithms/Algorithms.h index 7f27784b..47a0a92b 100644 --- a/Grid/algorithms/Algorithms.h +++ b/Grid/algorithms/Algorithms.h @@ -54,6 +54,7 @@ NAMESPACE_CHECK(BiCGSTAB); #include #include #include +#include #include #include #include diff --git a/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h b/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h index cd2e4374..31ac55e0 100644 --- a/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h +++ b/Grid/algorithms/iterative/ConjugateGradientMixedPrec.h @@ -49,6 +49,7 @@ NAMESPACE_BEGIN(Grid); Integer TotalInnerIterations; //Number of inner CG iterations Integer TotalOuterIterations; //Number of restarts Integer TotalFinalStepIterations; //Number of CG iterations in final patch-up step + RealD TrueResidual; //Option to speed up *inner single precision* solves using a LinearFunction that produces a guess LinearFunction *guesser; @@ -68,6 +69,7 @@ NAMESPACE_BEGIN(Grid); } void operator() (const FieldD &src_d_in, FieldD &sol_d){ + std::cout << GridLogMessage << "MixedPrecisionConjugateGradient: Starting mixed precision CG with outer tolerance " << Tolerance << " and inner tolerance " << InnerTolerance << std::endl; TotalInnerIterations = 0; GridStopWatch TotalTimer; @@ -97,6 +99,7 @@ NAMESPACE_BEGIN(Grid); FieldF sol_f(SinglePrecGrid); sol_f.Checkerboard() = cb; + std::cout< CG_f(inner_tol, MaxInnerIterations); CG_f.ErrorOnNoConverge = false; @@ -130,6 +133,7 @@ NAMESPACE_BEGIN(Grid); (*guesser)(src_f, sol_f); //Inner CG + std::cout< CG_d(Tolerance, MaxInnerIterations); CG_d(Linop_d, src_d_in, sol_d); TotalFinalStepIterations = CG_d.IterationsToComplete; + TrueResidual = CG_d.TrueResidual; TotalTimer.Stop(); std::cout< TrueResidualShift; - ConjugateGradientMultiShift(Integer maxit,MultiShiftFunction &_shifts) : + ConjugateGradientMultiShift(Integer maxit, const MultiShiftFunction &_shifts) : MaxIterations(maxit), shifts(_shifts) { @@ -182,6 +182,9 @@ public: for(int s=0;s +class ShiftedLinop: public LinearOperatorBase{ +public: + LinearOperatorBase &linop_base; + RealD shift; + + ShiftedLinop(LinearOperatorBase &_linop_base, RealD _shift): linop_base(_linop_base), shift(_shift){} + + void OpDiag (const Field &in, Field &out){ assert(0); } + void OpDir (const Field &in, Field &out,int dir,int disp){ assert(0); } + void OpDirAll (const Field &in, std::vector &out){ assert(0); } + + void Op (const Field &in, Field &out){ assert(0); } + void AdjOp (const Field &in, Field &out){ assert(0); } + + void HermOp(const Field &in, Field &out){ + linop_base.HermOp(in, out); + axpy(out, shift, in, out); + } + + void HermOpAndNorm(const Field &in, Field &out,RealD &n1,RealD &n2){ + HermOp(in,out); + ComplexD dot = innerProduct(in,out); + n1=real(dot); + n2=norm2(out); + } +}; +}; + + +template::value == 2, int>::type = 0, + typename std::enable_if< getPrecision::value == 1, int>::type = 0> +class ConjugateGradientMultiShiftMixedPrec : public OperatorMultiFunction, + public OperatorFunction +{ +public: + + using OperatorFunction::operator(); + + RealD Tolerance; + Integer MaxIterations; + Integer IterationsToComplete; //Number of iterations the CG took to finish. Filled in upon completion + std::vector IterationsToCompleteShift; // Iterations for this shift + int verbose; + MultiShiftFunction shifts; + std::vector TrueResidualShift; + + int ReliableUpdateFreq; //number of iterations between reliable updates + + GridBase* SinglePrecGrid; //Grid for single-precision fields + LinearOperatorBase &Linop_f; //single precision + + ConjugateGradientMultiShiftMixedPrec(Integer maxit, const MultiShiftFunction &_shifts, + GridBase* _SinglePrecGrid, LinearOperatorBase &_Linop_f, + int _ReliableUpdateFreq + ) : + MaxIterations(maxit), shifts(_shifts), SinglePrecGrid(_SinglePrecGrid), Linop_f(_Linop_f), ReliableUpdateFreq(_ReliableUpdateFreq) + { + verbose=1; + IterationsToCompleteShift.resize(_shifts.order); + TrueResidualShift.resize(_shifts.order); + } + + void operator() (LinearOperatorBase &Linop, const FieldD &src, FieldD &psi) + { + GridBase *grid = src.Grid(); + int nshift = shifts.order; + std::vector results(nshift,grid); + (*this)(Linop,src,results,psi); + } + void operator() (LinearOperatorBase &Linop, const FieldD &src, std::vector &results, FieldD &psi) + { + int nshift = shifts.order; + + (*this)(Linop,src,results); + + psi = shifts.norm*src; + for(int i=0;i &Linop_d, const FieldD &src_d, std::vector &psi_d) + { + GridBase *DoublePrecGrid = src_d.Grid(); + + //////////////////////////////////////////////////////////////////////// + // Convenience references to the info stored in "MultiShiftFunction" + //////////////////////////////////////////////////////////////////////// + int nshift = shifts.order; + + std::vector &mass(shifts.poles); // Make references to array in "shifts" + std::vector &mresidual(shifts.tolerances); + std::vector alpha(nshift,1.0); + + //Double precision search directions + FieldD p_d(DoublePrecGrid); + std::vector ps_d(nshift, DoublePrecGrid);// Search directions (double precision) + + FieldD tmp_d(DoublePrecGrid); + FieldD r_d(DoublePrecGrid); + FieldD mmp_d(DoublePrecGrid); + + assert(psi_d.size()==nshift); + assert(mass.size()==nshift); + assert(mresidual.size()==nshift); + + // dynamic sized arrays on stack; 2d is a pain with vector + RealD bs[nshift]; + RealD rsq[nshift]; + RealD z[nshift][2]; + int converged[nshift]; + + const int primary =0; + + //Primary shift fields CG iteration + RealD a,b,c,d; + RealD cp,bp,qq; //prev + + // Matrix mult fields + FieldF r_f(SinglePrecGrid); + FieldF p_f(SinglePrecGrid); + FieldF tmp_f(SinglePrecGrid); + FieldF mmp_f(SinglePrecGrid); + FieldF src_f(SinglePrecGrid); + precisionChange(src_f, src_d); + + // Check lightest mass + for(int s=0;s= mass[primary] ); + converged[s]=0; + } + + // Wire guess to zero + // Residuals "r" are src + // First search direction "p" is also src + cp = norm2(src_d); + + // Handle trivial case of zero src. + if( cp == 0. ){ + for(int s=0;s= rsq[s]){ + CleanupTimer.Start(); + std::cout< Linop_shift_d(Linop_d, mass[s]); + ConjugateGradientMultiShiftMixedPrecSupport::ShiftedLinop Linop_shift_f(Linop_f, mass[s]); + + MixedPrecisionConjugateGradient cg(mresidual[s], MaxIterations, MaxIterations, SinglePrecGrid, Linop_shift_f, Linop_shift_d); + cg(src_d, psi_d[s]); + + TrueResidualShift[s] = cg.TrueResidual; + CleanupTimer.Stop(); + } + } + + std::cout << GridLogMessage << "ConjugateGradientMultiShiftMixedPrec: Time Breakdown for body"< nbasis ) eresid = eresid*_coarse_relax_tol; + std::cout.precision(13); std::cout< nbasis ) eresid = eresid*_coarse_relax_tol; if( (vv on the coarse grid. This function orthnormalizes the fine-grid subspace + //vectors under the block inner product. This step must be performed after computing the fine grid + //eigenvectors and before computing the coarse grid eigenvectors. void Orthogonalise(void ) { CoarseScalar InnerProd(_CoarseGrid); std::cout << GridLogMessage <<" Gramm-Schmidt pass 1"< Cheby(cheby_op); - ProjectedHermOp Op(_FineOp,subspace); - ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,subspace); + Chebyshev Cheby(cheby_op); //Chebyshev of fine operator on fine grid + ProjectedHermOp Op(_FineOp,subspace); //Fine operator on coarse grid with intermediate fine grid conversion + ProjectedFunctionHermOp ChebyOp (Cheby,_FineOp,subspace); //Chebyshev of fine operator on coarse grid with intermediate fine grid conversion ////////////////////////////////////////////////////////////////////////////////////////////////// // create a smoother and see if we can get a cheap convergence test and smooth inside the IRL ////////////////////////////////////////////////////////////////////////////////////////////////// - Chebyshev ChebySmooth(cheby_smooth); - ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); + Chebyshev ChebySmooth(cheby_smooth); //lower order Chebyshev of fine operator on fine grid used to smooth regenerated eigenvectors + ImplicitlyRestartedLanczosSmoothedTester ChebySmoothTester(ChebyOp,ChebySmooth,_FineOp,subspace,relax); evals_coarse.resize(Nm); evec_coarse.resize(Nm,_CoarseGrid); CoarseField src(_CoarseGrid); src=1.0; + //Note the "tester" here is also responsible for generating the fine grid eigenvalues which are output into the "evals_coarse" array ImplicitlyRestartedLanczos IRL(ChebyOp,ChebyOp,ChebySmoothTester,Nstop,Nk,Nm,resid,MaxIt,betastp,MinRes); int Nconv=0; IRL.calc(evals_coarse,evec_coarse,src,Nconv,false); @@ -405,6 +427,14 @@ public: std::cout << i << " Coarse eval = " << evals_coarse[i] << std::endl; } } + + //Get the fine eigenvector 'i' by reconstruction + void getFineEvecEval(FineField &evec, RealD &eval, const int i) const{ + blockPromote(evec_coarse[i],evec,subspace); + eval = evals_coarse[i]; + } + + }; NAMESPACE_END(Grid); diff --git a/Grid/algorithms/iterative/PowerMethod.h b/Grid/algorithms/iterative/PowerMethod.h index 6aa8e923..027ea68c 100644 --- a/Grid/algorithms/iterative/PowerMethod.h +++ b/Grid/algorithms/iterative/PowerMethod.h @@ -29,6 +29,8 @@ template class PowerMethod RealD vnum = real(innerProduct(src_n,tmp)); // HermOp. RealD vden = norm2(src_n); RealD na = vnum/vden; + + std::cout << GridLogIterative << "PowerMethod: Current approximation of largest eigenvalue " << na << std::endl; if ( (fabs(evalMaxApprox/na - 1.0) < 0.001) || (i==_MAX_ITER_EST_-1) ) { evalMaxApprox = na; diff --git a/Grid/parallelIO/NerscIO.h b/Grid/parallelIO/NerscIO.h index 99011e25..88278131 100644 --- a/Grid/parallelIO/NerscIO.h +++ b/Grid/parallelIO/NerscIO.h @@ -39,9 +39,11 @@ using namespace Grid; //////////////////////////////////////////////////////////////////////////////// class NerscIO : public BinaryIO { public: - typedef Lattice GaugeField; + // Enable/disable exiting if the plaquette in the header does not match the value computed (default true) + static bool & exitOnReadPlaquetteMismatch(){ static bool v=true; return v; } + static inline void truncate(std::string file){ std::ofstream fout(file,std::ios::out); } @@ -198,7 +200,7 @@ public: std::cerr << " nersc_csum " < class GparityWilsonImpl : public ConjugateGaugeImpl > { public: @@ -113,7 +125,7 @@ public: || ((distance== 1)&&(icoor[direction]==1)) || ((distance==-1)&&(icoor[direction]==0)); - permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu]; //only if we are going around the world + permute_lane = permute_lane && SE->_around_the_world && St.parameters.twists[mmu] && mmu < Nd-1; //only if we are going around the world in a spatial direction //Apply the links int f_upper = permute_lane ? 1 : 0; @@ -139,10 +151,10 @@ public: assert((distance == 1) || (distance == -1)); // nearest neighbour stencil hard code assert((sl == 1) || (sl == 2)); - if ( SE->_around_the_world && St.parameters.twists[mmu] ) { - + //If this site is an global boundary site, perform the G-parity flavor twist + if ( mmu < Nd-1 && SE->_around_the_world && St.parameters.twists[mmu] ) { if ( sl == 2 ) { - + //Only do the twist for lanes on the edge of the physical node ExtractBuffer vals(Nsimd); extract(chi,vals); @@ -197,6 +209,19 @@ public: reg = memory; } + + //Poke 'poke_f0' onto flavor 0 and 'poke_f1' onto flavor 1 in direction mu of the doubled gauge field Uds + inline void pokeGparityDoubledGaugeField(DoubledGaugeField &Uds, const GaugeLinkField &poke_f0, const GaugeLinkField &poke_f1, const int mu){ + autoView(poke_f0_v, poke_f0, CpuRead); + autoView(poke_f1_v, poke_f1, CpuRead); + autoView(Uds_v, Uds, CpuWrite); + thread_foreach(ss,poke_f0_v,{ + Uds_v[ss](0)(mu) = poke_f0_v[ss](); + Uds_v[ss](1)(mu) = poke_f1_v[ss](); + }); + } + + inline void DoubleStore(GridBase *GaugeGrid,DoubledGaugeField &Uds,const GaugeField &Umu) { conformable(Uds.Grid(),GaugeGrid); @@ -207,14 +232,19 @@ public: GaugeLinkField Uconj(GaugeGrid); Lattice > coor(GaugeGrid); - - for(int mu=0;mu(Umu,mu); Uconj = conjugate(U); + // Implement the isospin rotation sign on the boundary between f=1 and f=0 // This phase could come from a simple bc 1,1,-1,1 .. int neglink = GaugeGrid->GlobalDimensions()[mu]-1; if ( Params.twists[mu] ) { @@ -229,7 +259,7 @@ public: thread_foreach(ss,U_v,{ Uds_v[ss](0)(mu) = U_v[ss](); Uds_v[ss](1)(mu) = Uconj_v[ss](); - }); + }); } U = adj(Cshift(U ,mu,-1)); // correct except for spanning the boundary @@ -260,6 +290,38 @@ public: }); } } + + { //periodic / antiperiodic temporal BCs + int mu = Nd-1; + int L = GaugeGrid->GlobalDimensions()[mu]; + int Lmu = L - 1; + + LatticeCoordinate(coor, mu); + + U = PeekIndex(Umu, mu); //Get t-directed links + + GaugeLinkField *Upoke = &U; + + if(Params.twists[mu]){ //antiperiodic + Utmp = where(coor == Lmu, -U, U); + Upoke = &Utmp; + } + + Uconj = conjugate(*Upoke); //second flavor interacts with conjugate links + pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu); + + //Get the barrel-shifted field + Utmp = adj(Cshift(U, mu, -1)); //is a forward shift! + Upoke = &Utmp; + + if(Params.twists[mu]){ + U = where(coor == 0, -Utmp, Utmp); //boundary phase + Upoke = &U; + } + + Uconj = conjugate(*Upoke); + pokeGparityDoubledGaugeField(Uds, *Upoke, Uconj, mu + 4); + } } inline void InsertForce4D(GaugeField &mat, FermionField &Btilde, FermionField &A, int mu) { @@ -298,28 +360,48 @@ public: inline void extractLinkField(std::vector &mat, DoubledGaugeField &Uds){ assert(0); } - + inline void InsertForce5D(GaugeField &mat, FermionField &Btilde, FermionField Ã, int mu) { - - int Ls = Btilde.Grid()->_fdimensions[0]; - - GaugeLinkField tmp(mat.Grid()); - tmp = Zero(); + int Ls=Btilde.Grid()->_fdimensions[0]; + { - autoView( tmp_v , tmp, CpuWrite); - autoView( Atilde_v , Atilde, CpuRead); - autoView( Btilde_v , Btilde, CpuRead); - thread_for(ss,tmp.Grid()->oSites(),{ - for (int s = 0; s < Ls; s++) { - int sF = s + Ls * ss; - auto ttmp = traceIndex(outerProduct(Btilde_v[sF], Atilde_v[sF])); - tmp_v[ss]() = tmp_v[ss]() + ttmp(0, 0) + conjugate(ttmp(1, 1)); - } - }); + GridBase *GaugeGrid = mat.Grid(); + Lattice > coor(GaugeGrid); + + if( Params.twists[mu] ){ + LatticeCoordinate(coor,mu); + } + + autoView( mat_v , mat, AcceleratorWrite); + autoView( Btilde_v , Btilde, AcceleratorRead); + autoView( Atilde_v , Atilde, AcceleratorRead); + accelerator_for(sss,mat.Grid()->oSites(), FermionField::vector_type::Nsimd(),{ + int sU=sss; + typedef decltype(coalescedRead(mat_v[sU](mu)() )) ColorMatrixType; + ColorMatrixType sum; + zeroit(sum); + for(int s=0;s(mat, tmp, mu); - return; } + + + + }; diff --git a/tests/core/Test_gparity.cc b/tests/core/Test_gparity.cc index b2068901..5bf98ba6 100644 --- a/tests/core/Test_gparity.cc +++ b/tests/core/Test_gparity.cc @@ -55,13 +55,17 @@ static_assert(same_vComplex == 1, "Dirac Operators must have same underlying SIM int main (int argc, char ** argv) { int nu = 0; - + int tbc_aprd = 0; //use antiperiodic BCs in the time direction? + Grid_init(&argc,&argv); for(int i=1;i> nu; std::cout << GridLogMessage << "Set Gparity direction to " << nu << std::endl; + }else if(std::string(argv[i]) == "--Tbc-APRD"){ + tbc_aprd = 1; + std::cout << GridLogMessage << "Using antiperiodic BCs in the time direction" << std::endl; } } @@ -155,13 +159,18 @@ int main (int argc, char ** argv) //Coordinate grid for reference LatticeInteger xcoor_1f5(FGrid_1f); - LatticeCoordinate(xcoor_1f5,1+nu); + LatticeCoordinate(xcoor_1f5,1+nu); //note '1+nu'! This is because for 5D fields the s-direction is direction 0 Replicate(src,src_1f); src_1f = where( xcoor_1f5 >= Integer(L), 2.0*src_1f,src_1f ); RealD mass=0.0; RealD M5=1.8; - StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS); + + //Standard Dirac op + AcceleratorVector bc_std(Nd, 1.0); + if(tbc_aprd) bc_std[Nd-1] = -1.; //antiperiodic time BC + StandardDiracOp::ImplParams std_params(bc_std); + StandardDiracOp Ddwf(Umu_1f,*FGrid_1f,*FrbGrid_1f,*UGrid_1f,*UrbGrid_1f,mass,M5 DOP_PARAMS, std_params); StandardFermionField src_o_1f(FrbGrid_1f); StandardFermionField result_o_1f(FrbGrid_1f); @@ -172,9 +181,11 @@ int main (int argc, char ** argv) ConjugateGradient CG(1.0e-8,10000); CG(HermOpEO,src_o_1f,result_o_1f); - // const int nu = 3; + //Gparity Dirac op std::vector twists(Nd,0); twists[nu] = 1; + if(tbc_aprd) twists[Nd-1] = 1; + GparityDiracOp::ImplParams params; params.twists = twists; GparityDiracOp GPDdwf(Umu_2f,*FGrid_2f,*FrbGrid_2f,*UGrid_2f,*UrbGrid_2f,mass,M5 DOP_PARAMS,params); @@ -271,8 +282,11 @@ int main (int argc, char ** argv) std::cout << "2f cb "<(result_o_2f,0); - res1o = PeekIndex<0>(result_o_2f,1); + res0o = PeekIndex<0>(result_o_2f,0); //flavor 0, odd cb + res1o = PeekIndex<0>(result_o_2f,1); //flavor 1, odd cb std::cout << "res cb "<= Integer(L), replica1,replica0 ); replica0 = Zero(); setCheckerboard(replica0,result_o_1f); - std::cout << "Norm2 solutions is " < twists(Nd,0); // twists[nu] = 1; - // GparityDomainWallFermionR::ImplParams params; params.twists = twists; - // GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); - // DomainWallFermionR Dw (U, Grid,RBGrid,mass,M5); - - const int nu = 3; + const int nu = 0; //gparity direction std::vector twists(Nd,0); twists[nu] = 1; + twists[Nd-1] = 1; //antiperiodic in time GparityDomainWallFermionR::ImplParams params; params.twists = twists; - - /* - params.boundary_phases[0] = 1.0; - params.boundary_phases[1] = 1.0; - params.boundary_phases[2] = 1.0; - params.boundary_phases[3] =- 1.0; - */ - + GparityDomainWallFermionR Dw(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); Dw.M (phi,Mphi); diff --git a/tests/forces/Test_gpdwf_force.cc b/tests/forces/Test_gpdwf_force.cc index d6744080..af1ce82b 100644 --- a/tests/forces/Test_gpdwf_force.cc +++ b/tests/forces/Test_gpdwf_force.cc @@ -71,8 +71,10 @@ int main (int argc, char ** argv) RealD mass=0.01; RealD M5=1.8; - const int nu = 3; - std::vector twists(Nd,0); twists[nu] = 1; + const int nu = 1; + std::vector twists(Nd,0); + twists[nu] = 1; + twists[3] = 1; GparityDomainWallFermionR::ImplParams params; params.twists = twists; GparityDomainWallFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,params); Ddwf.M (phi,Mphi); diff --git a/tests/forces/Test_gpwilson_force.cc b/tests/forces/Test_gpwilson_force.cc index d731f27a..7ab2ddeb 100644 --- a/tests/forces/Test_gpwilson_force.cc +++ b/tests/forces/Test_gpwilson_force.cc @@ -64,8 +64,12 @@ int main (int argc, char ** argv) //////////////////////////////////// RealD mass=0.01; - const int nu = 3; - std::vector twists(Nd,0); twists[nu] = 1; + const int nu = 1; + const int Lnu=latt_size[nu]; + + std::vector twists(Nd,0); + twists[nu] = 1; + twists[3]=1; GparityWilsonFermionR::ImplParams params; params.twists = twists; GparityWilsonFermionR Wil(U,*UGrid,*UrbGrid,mass,params); Wil.M (phi,Mphi); diff --git a/tests/lanczos/Test_dwf_lanczos.cc b/tests/lanczos/Test_dwf_lanczos.cc index 00d29ec0..1fe29bb2 100644 --- a/tests/lanczos/Test_dwf_lanczos.cc +++ b/tests/lanczos/Test_dwf_lanczos.cc @@ -31,14 +31,38 @@ using namespace std; using namespace Grid; ; -typedef typename GparityDomainWallFermionR::FermionField FermionField; +template +struct Setup{}; -RealD AllZero(RealD x){ return 0.;} +template<> +struct Setup{ + static GparityMobiusFermionR* getAction(LatticeGaugeField &Umu, + GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){ + RealD mass=0.01; + RealD M5=1.8; + RealD mob_b=1.5; + GparityMobiusFermionD ::ImplParams params; + std::vector twists({1,1,1,0}); + params.twists = twists; + return new GparityMobiusFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); + } +}; -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); +template<> +struct Setup{ + static DomainWallFermionR* getAction(LatticeGaugeField &Umu, + GridCartesian* FGrid, GridRedBlackCartesian* FrbGrid, GridCartesian* UGrid, GridRedBlackCartesian* UrbGrid){ + RealD mass=0.01; + RealD M5=1.8; + return new DomainWallFermionR(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); + } +}; + + +template +void run(){ + typedef typename Action::FermionField FermionField; const int Ls=8; GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi()); @@ -56,24 +80,10 @@ int main (int argc, char ** argv) LatticeGaugeField Umu(UGrid); SU::HotConfiguration(RNG4, Umu); - std::vector U(4,UGrid); - for(int mu=0;mu(Umu,mu); - } - - RealD mass=0.01; - RealD M5=1.8; - RealD mob_b=1.5; -// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5); - GparityMobiusFermionD ::ImplParams params; - std::vector twists({1,1,1,0}); - params.twists = twists; - GparityMobiusFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,mob_b,mob_b-1.,params); - -// MdagMLinearOperator HermOp(Ddwf); -// SchurDiagTwoOperator HermOp(Ddwf); - SchurDiagTwoOperator HermOp(Ddwf); -// SchurDiagMooeeOperator HermOp(Ddwf); + Action *action = Setup::getAction(Umu,FGrid,FrbGrid,UGrid,UrbGrid); + + //MdagMLinearOperator HermOp(Ddwf); + SchurDiagTwoOperator HermOp(*action); const int Nstop = 30; const int Nk = 40; @@ -90,8 +100,7 @@ int main (int argc, char ** argv) PlainHermOp Op (HermOp); ImplicitlyRestartedLanczos IRL(OpCheby,Op,Nstop,Nk,Nm,resid,MaxIt); - - + std::vector eval(Nm); FermionField src(FrbGrid); gaussian(RNG5rb,src); @@ -103,6 +112,28 @@ int main (int argc, char ** argv) int Nconv; IRL.calc(eval,evec,src,Nconv); + delete action; +} + +int main (int argc, char ** argv) +{ + Grid_init(&argc,&argv); + std::string action = "GparityMobius"; + for(int i=1;i(); + }else if(action == "DWF"){ + run(); + }else{ + std::cout << "Unknown action" << std::endl; + exit(1); + } + Grid_finalize(); } diff --git a/tests/solver/Test_dwf_multishift_mixedprec.cc b/tests/solver/Test_dwf_multishift_mixedprec.cc new file mode 100644 index 00000000..bdede459 --- /dev/null +++ b/tests/solver/Test_dwf_multishift_mixedprec.cc @@ -0,0 +1,184 @@ + /************************************************************************************* + + Grid physics library, www.github.com/paboyle/Grid + + Source file: ./tests/Test_dwf_multishift_mixedprec.cc + + Copyright (C) 2015 + +Author: Christopher Kelly + + 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 + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + This program is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License along + with this program; if not, write to the Free Software Foundation, Inc., + 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + + See the full license in the file "LICENSE" in the top level distribution directory + *************************************************************************************/ + /* END LEGAL */ +#include + +using namespace Grid; + +template +void run_test(int argc, char ** argv, const typename SpeciesD::ImplParams ¶ms){ + const int Ls = 16; + GridCartesian* UGrid_d = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexD::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid_d = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_d); + GridCartesian* FGrid_d = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_d); + GridRedBlackCartesian* FrbGrid_d = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_d); + + GridCartesian* UGrid_f = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd, vComplexF::Nsimd()), GridDefaultMpi()); + GridRedBlackCartesian* UrbGrid_f = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid_f); + GridCartesian* FGrid_f = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid_f); + GridRedBlackCartesian* FrbGrid_f = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid_f); + + typedef typename SpeciesD::FermionField FermionFieldD; + typedef typename SpeciesF::FermionField FermionFieldF; + + std::vector seeds4({1, 2, 3, 4}); + std::vector seeds5({5, 6, 7, 8}); + GridParallelRNG RNG5(FGrid_d); + RNG5.SeedFixedIntegers(seeds5); + GridParallelRNG RNG4(UGrid_d); + RNG4.SeedFixedIntegers(seeds4); + + FermionFieldD src_d(FGrid_d); + random(RNG5, src_d); + + LatticeGaugeFieldD Umu_d(UGrid_d); + + //CPS-created G-parity ensembles have a factor of 2 error in the plaquette that causes the read to fail unless we workaround it + bool gparity_plaquette_fix = false; + for(int i=1;i(Umu_d, metadata, file); + + if(gparity_plaquette_fix){ + metadata.plaquette *= 2.; //correct header value + + //Get the true plaquette + FieldMetaData tmp; + GaugeStatisticsType gs; gs(Umu_d, tmp); + + std::cout << "After correction: plaqs " << tmp.plaquette << " " << metadata.plaquette << std::endl; + assert(fabs(tmp.plaquette -metadata.plaquette ) < 1.0e-5 ); + } + + cfg_loaded=true; + break; + } + } + + if(!cfg_loaded) + SU::HotConfiguration(RNG4, Umu_d); + + LatticeGaugeFieldF Umu_f(UGrid_f); + precisionChange(Umu_f, Umu_d); + + std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt() << " Ls: " << Ls << std::endl; + + RealD mass = 0.01; + RealD M5 = 1.8; + SpeciesD Ddwf_d(Umu_d, *FGrid_d, *FrbGrid_d, *UGrid_d, *UrbGrid_d, mass, M5, params); + SpeciesF Ddwf_f(Umu_f, *FGrid_f, *FrbGrid_f, *UGrid_f, *UrbGrid_f, mass, M5, params); + + FermionFieldD src_o_d(FrbGrid_d); + pickCheckerboard(Odd, src_o_d, src_d); + + SchurDiagMooeeOperator HermOpEO_d(Ddwf_d); + SchurDiagMooeeOperator HermOpEO_f(Ddwf_f); + + AlgRemez remez(1e-4, 64, 50); + int order = 15; + remez.generateApprox(order, 1, 2); //sqrt + + MultiShiftFunction shifts(remez, 1e-10, false); + + int relup_freq = 50; + double t1=usecond(); + ConjugateGradientMultiShiftMixedPrec mcg(10000, shifts, FrbGrid_f, HermOpEO_f, relup_freq); + + std::vector results_o_d(order, FrbGrid_d); + mcg(HermOpEO_d, src_o_d, results_o_d); + double t2=usecond(); + + //Crosscheck double and mixed prec results + ConjugateGradientMultiShift dmcg(10000, shifts); + std::vector results_o_d_2(order, FrbGrid_d); + dmcg(HermOpEO_d, src_o_d, results_o_d_2); + double t3=usecond(); + + std::cout << GridLogMessage << "Comparison of mixed prec results to double prec results |mixed - double|^2 :" << std::endl; + FermionFieldD tmp(FrbGrid_d); + for(int i=0;i= 0 && gpdir <= 2); //spatial! + gparity = true; + } + } + if(gparity){ + std::cout << "Running test with G-parity BCs in " << gpdir << " direction" << std::endl; + GparityWilsonImplParams params; + params.twists[gpdir] = 1; + + std::vector conj_dirs(Nd,0); + conj_dirs[gpdir] = 1; + ConjugateGimplD::setDirections(conj_dirs); + + run_test(argc,argv,params); + }else{ + std::cout << "Running test with periodic BCs" << std::endl; + WilsonImplParams params; + run_test(argc,argv,params); + } + + Grid_finalize(); +} From 1ad54d049d2dc24575714c6e02f7ef6a01b3ddaa Mon Sep 17 00:00:00 2001 From: Christopher Kelly Date: Thu, 2 Jun 2022 15:30:41 -0400 Subject: [PATCH 4/4] To PeriodicBC and ConjugateBC, added a new function "CshiftLink" which performs a boundary-aware C-shift of links or products of links. For the latter, the links crossing the global boundary are complex-conjugated. To the gauge implementations, added CshiftLink functions calling into the appropriate operation for the BC in a given direction. GaugeTransform, FourierAcceleratedGaugeFixer and WilsonLoops::FieldStrength no longer implicitly assume periodic boundary conditions; instead the shifted link is obtained using CshiftLink and is aware of the gauge implementation. Added an assert-check to ensure that the gauge fixing converges within the specified number of steps. Added functionality to compute the timeslice averaged plaquette Added functionality to compute the 5LI topological charge and timeslice topological charge Added a check of the properties of the charge conjugation matrix C=-gamma_2 gamma_4 to Test_gamma Fixed const correctness for Replicate Modified Test_fft_gfix to support either conjugate or periodic BCs, optionally disabling Fourier-accelerated gauge fixing, and tuning of alpha using cmdline options --- Grid/lattice/Lattice_transfer.h | 2 +- Grid/qcd/action/gauge/GaugeImplementations.h | 38 +++ Grid/qcd/utils/CovariantCshift.h | 41 +++ Grid/qcd/utils/GaugeFix.h | 62 +++-- Grid/qcd/utils/SUn.h | 24 +- Grid/qcd/utils/WilsonLoops.h | 251 ++++++++++++++++++- examples/Example_Laplacian.cc | 2 +- tests/core/Test_fft_gfix.cc | 187 ++++++++++---- tests/core/Test_gamma.cc | 60 +++++ 9 files changed, 572 insertions(+), 95 deletions(-) diff --git a/Grid/lattice/Lattice_transfer.h b/Grid/lattice/Lattice_transfer.h index ef489ea6..aee55e93 100644 --- a/Grid/lattice/Lattice_transfer.h +++ b/Grid/lattice/Lattice_transfer.h @@ -855,7 +855,7 @@ void ExtractSliceLocal(Lattice &lowDim,const Lattice & higherDim,int template -void Replicate(Lattice &coarse,Lattice & fine) +void Replicate(const Lattice &coarse,Lattice & fine) { typedef typename vobj::scalar_object sobj; diff --git a/Grid/qcd/action/gauge/GaugeImplementations.h b/Grid/qcd/action/gauge/GaugeImplementations.h index 16147c77..f518b236 100644 --- a/Grid/qcd/action/gauge/GaugeImplementations.h +++ b/Grid/qcd/action/gauge/GaugeImplementations.h @@ -69,6 +69,11 @@ public: return PeriodicBC::ShiftStaple(Link,mu); } + //Same as Cshift for periodic BCs + static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){ + return PeriodicBC::CshiftLink(Link,mu,shift); + } + static inline bool isPeriodicGaugeField(void) { return true; } }; @@ -110,6 +115,11 @@ public: return PeriodicBC::CovShiftBackward(Link, mu, field); } + //If mu is a conjugate BC direction + //Out(x) = U^dag_\mu(x-mu) | x_\mu != 0 + // = U^T_\mu(L-1) | x_\mu == 0 + //else + //Out(x) = U^dag_\mu(x-mu mod L) static inline GaugeLinkField CovShiftIdentityBackward(const GaugeLinkField &Link, int mu) { @@ -129,6 +139,13 @@ public: return PeriodicBC::CovShiftIdentityForward(Link,mu); } + + //If mu is a conjugate BC direction + //Out(x) = S_\mu(x+mu) | x_\mu != L-1 + // = S*_\mu(x+mu) | x_\mu == L-1 + //else + //Out(x) = S_\mu(x+mu mod L) + //Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices static inline GaugeLinkField ShiftStaple(const GaugeLinkField &Link, int mu) { assert(_conjDirs.size() == Nd); @@ -138,6 +155,27 @@ public: return PeriodicBC::ShiftStaple(Link,mu); } + //Boundary-aware C-shift of gauge links / gauge transformation matrices + //For conjugate BC direction + //shift = 1 + //Out(x) = U_\mu(x+\hat\mu) | x_\mu != L-1 + // = U*_\mu(0) | x_\mu == L-1 + //shift = -1 + //Out(x) = U_\mu(x-mu) | x_\mu != 0 + // = U*_\mu(L-1) | x_\mu == 0 + //else + //shift = 1 + //Out(x) = U_\mu(x+\hat\mu mod L) + //shift = -1 + //Out(x) = U_\mu(x-\hat\mu mod L) + static inline GaugeLinkField CshiftLink(const GaugeLinkField &Link, int mu, int shift){ + assert(_conjDirs.size() == Nd); + if(_conjDirs[mu]) + return ConjugateBC::CshiftLink(Link,mu,shift); + else + return PeriodicBC::CshiftLink(Link,mu,shift); + } + static inline void setDirections(std::vector &conjDirs) { _conjDirs=conjDirs; } static inline std::vector getDirections(void) { return _conjDirs; } static inline bool isPeriodicGaugeField(void) { return false; } diff --git a/Grid/qcd/utils/CovariantCshift.h b/Grid/qcd/utils/CovariantCshift.h index 6c70706f..79cf8e0f 100644 --- a/Grid/qcd/utils/CovariantCshift.h +++ b/Grid/qcd/utils/CovariantCshift.h @@ -88,6 +88,12 @@ namespace PeriodicBC { return CovShiftBackward(Link,mu,arg); } + //Boundary-aware C-shift of gauge links / gauge transformation matrices + template Lattice + CshiftLink(const Lattice &Link, int mu, int shift) + { + return Cshift(Link, mu, shift); + } } @@ -158,6 +164,9 @@ namespace ConjugateBC { // std::cout<<"Gparity::CovCshiftBackward mu="< Lattice CovShiftIdentityBackward(const Lattice &Link, int mu) { GridBase *grid = Link.Grid(); @@ -176,6 +185,9 @@ namespace ConjugateBC { return Link; } + //Out(x) = S_\mu(x+\hat\mu) | x_\mu != L-1 + // = S*_\mu(0) | x_\mu == L-1 + //Note: While this is used for Staples it is also applicable for shifting gauge links or gauge transformation matrices template Lattice ShiftStaple(const Lattice &Link, int mu) { @@ -208,6 +220,35 @@ namespace ConjugateBC { return CovShiftBackward(Link,mu,arg); } + //Boundary-aware C-shift of gauge links / gauge transformation matrices + //shift = 1 + //Out(x) = U_\mu(x+\hat\mu) | x_\mu != L-1 + // = U*_\mu(0) | x_\mu == L-1 + //shift = -1 + //Out(x) = U_\mu(x-mu) | x_\mu != 0 + // = U*_\mu(L-1) | x_\mu == 0 + template Lattice + CshiftLink(const Lattice &Link, int mu, int shift) + { + GridBase *grid = Link.Grid(); + int Lmu = grid->GlobalDimensions()[mu] - 1; + + Lattice> coor(grid); + LatticeCoordinate(coor, mu); + + Lattice tmp(grid); + if(shift == 1){ + tmp = Cshift(Link, mu, 1); + tmp = where(coor == Lmu, conjugate(tmp), tmp); + return tmp; + }else if(shift == -1){ + tmp = Link; + tmp = where(coor == Lmu, conjugate(tmp), tmp); + return Cshift(tmp, mu, -1); + }else assert(0 && "Invalid shift value"); + return tmp; //shuts up the compiler fussing about the return type + } + } diff --git a/Grid/qcd/utils/GaugeFix.h b/Grid/qcd/utils/GaugeFix.h index 2b3384da..c0bc2c83 100644 --- a/Grid/qcd/utils/GaugeFix.h +++ b/Grid/qcd/utils/GaugeFix.h @@ -40,27 +40,46 @@ public: typedef typename Gimpl::GaugeLinkField GaugeMat; typedef typename Gimpl::GaugeField GaugeLorentz; - static void GaugeLinkToLieAlgebraField(const std::vector &U,std::vector &A) { - for(int mu=0;mu &A,GaugeMat &dmuAmu,int orthog) { + + //The derivative of the Lie algebra field + static void DmuAmu(const std::vector &U, GaugeMat &dmuAmu,int orthog) { + GridBase* grid = U[0].Grid(); + GaugeMat Ax(grid); + GaugeMat Axm1(grid); + GaugeMat Utmp(grid); + dmuAmu=Zero(); for(int mu=0;mu &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) { + static Real SteepestDescentStep(std::vector &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) { GridBase *grid = U[0].Grid(); - std::vector A(Nd,grid); GaugeMat g(grid); - - GaugeLinkToLieAlgebraField(U,A); - ExpiAlphaDmuAmu(A,g,alpha,dmuAmu,orthog); - + ExpiAlphaDmuAmu(U,g,alpha,dmuAmu,orthog); Real vol = grid->gSites(); Real trG = TensorRemove(sum(trace(g))).real()/vol/Nc; xform = g*xform ; - SU::GaugeTransform(U,g); + SU::GaugeTransform(U,g); return trG; } - static Real FourierAccelSteepestDescentStep(std::vector &U,GaugeMat &xform,Real & alpha, GaugeMat & dmuAmu,int orthog) { + static Real FourierAccelSteepestDescentStep(std::vector &U,GaugeMat &xform, Real alpha, GaugeMat & dmuAmu,int orthog) { GridBase *grid = U[0].Grid(); @@ -157,11 +173,7 @@ public: GaugeMat g(grid); GaugeMat dmuAmu_p(grid); - std::vector A(Nd,grid); - - GaugeLinkToLieAlgebraField(U,A); - - DmuAmu(A,dmuAmu,orthog); + DmuAmu(U,dmuAmu,orthog); std::vector mask(Nd,1); for(int mu=0;mu::GaugeTransform(U,g); + SU::GaugeTransform(U,g); return trG; } - static void ExpiAlphaDmuAmu(const std::vector &A,GaugeMat &g,Real & alpha, GaugeMat &dmuAmu,int orthog) { + static void ExpiAlphaDmuAmu(const std::vector &U,GaugeMat &g, Real alpha, GaugeMat &dmuAmu,int orthog) { GridBase *grid = g.Grid(); Complex cialpha(0.0,-alpha); GaugeMat ciadmam(grid); - DmuAmu(A,dmuAmu,orthog); + DmuAmu(U,dmuAmu,orthog); ciadmam = dmuAmu*cialpha; SU::taExp(ciadmam,g); } diff --git a/Grid/qcd/utils/SUn.h b/Grid/qcd/utils/SUn.h index 675493b3..b9660c65 100644 --- a/Grid/qcd/utils/SUn.h +++ b/Grid/qcd/utils/SUn.h @@ -694,32 +694,32 @@ public: * Adjoint rep gauge xform */ - template - static void GaugeTransform( GaugeField &Umu, GaugeMat &g){ + template + static void GaugeTransform(typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){ GridBase *grid = Umu.Grid(); conformable(grid,g.Grid()); - GaugeMat U(grid); - GaugeMat ag(grid); ag = adj(g); + typename Gimpl::GaugeLinkField U(grid); + typename Gimpl::GaugeLinkField ag(grid); ag = adj(g); for(int mu=0;mu(Umu,mu); - U = g*U*Cshift(ag, mu, 1); + U = g*U*Gimpl::CshiftLink(ag, mu, 1); //BC-aware PokeIndex(Umu,U,mu); } } - template - static void GaugeTransform( std::vector &U, GaugeMat &g){ + template + static void GaugeTransform( std::vector &U, typename Gimpl::GaugeLinkField &g){ GridBase *grid = g.Grid(); - GaugeMat ag(grid); ag = adj(g); + typename Gimpl::GaugeLinkField ag(grid); ag = adj(g); for(int mu=0;mu - static void RandomGaugeTransform(GridParallelRNG &pRNG, GaugeField &Umu, GaugeMat &g){ + template + static void RandomGaugeTransform(GridParallelRNG &pRNG, typename Gimpl::GaugeField &Umu, typename Gimpl::GaugeLinkField &g){ LieRandomize(pRNG,g,1.0); - GaugeTransform(Umu,g); + GaugeTransform(Umu,g); } // Projects the algebra components a lattice matrix (of dimension ncol*ncol -1 ) diff --git a/Grid/qcd/utils/WilsonLoops.h b/Grid/qcd/utils/WilsonLoops.h index 0367c9fa..da1f5ac8 100644 --- a/Grid/qcd/utils/WilsonLoops.h +++ b/Grid/qcd/utils/WilsonLoops.h @@ -125,6 +125,56 @@ public: return sumplaq / vol / faces / Nc; // Nd , Nc dependent... FIXME } + ////////////////////////////////////////////////// + // sum over all spatial planes of plaquette + ////////////////////////////////////////////////// + static void siteSpatialPlaquette(ComplexField &Plaq, + const std::vector &U) { + ComplexField sitePlaq(U[0].Grid()); + Plaq = Zero(); + for (int mu = 1; mu < Nd-1; mu++) { + for (int nu = 0; nu < mu; nu++) { + traceDirPlaquette(sitePlaq, U, mu, nu); + Plaq = Plaq + sitePlaq; + } + } + } + + //////////////////////////////////// + // sum over all x,y,z and over all spatial planes of plaquette + ////////////////////////////////////////////////// + static std::vector timesliceSumSpatialPlaquette(const GaugeLorentz &Umu) { + std::vector U(Nd, Umu.Grid()); + // inefficient here + for (int mu = 0; mu < Nd; mu++) { + U[mu] = PeekIndex(Umu, mu); + } + + ComplexField Plaq(Umu.Grid()); + + siteSpatialPlaquette(Plaq, U); + typedef typename ComplexField::scalar_object sobj; + std::vector Tq; + sliceSum(Plaq, Tq, Nd-1); + + std::vector out(Tq.size()); + for(int t=0;t timesliceAvgSpatialPlaquette(const GaugeLorentz &Umu) { + std::vector sumplaq = timesliceSumSpatialPlaquette(Umu); + int Lt = Umu.Grid()->FullDimensions()[Nd-1]; + assert(sumplaq.size() == Lt); + double vol = Umu.Grid()->gSites() / Lt; + double faces = (1.0 * (Nd - 1)* (Nd - 2)) / 2.0; + for(int t=0;t(Umu, mu); // some redundant copies GaugeMat vu = v*u; //FS = 0.25*Ta(u*v + Cshift(vu, mu, -1)); - FS = (u*v + Cshift(vu, mu, -1)); + FS = (u*v + Gimpl::CshiftLink(vu, mu, -1)); FS = 0.125*(FS - adj(FS)); } - static Real TopologicalCharge(GaugeLorentz &U){ + static Real TopologicalCharge(const GaugeLorentz &U){ // 4d topological charge assert(Nd==4); // Bx = -iF(y,z), By = -iF(z,y), Bz = -iF(x,y) @@ -390,6 +440,203 @@ public: } + //Clover-leaf Wilson loop combination for arbitrary mu-extent M and nu extent N, mu >= nu + //cf https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 7 for 1x2 Wilson loop + //Clockwise ordering + static void CloverleafMxN(GaugeMat &FS, const GaugeMat &Umu, const GaugeMat &Unu, int mu, int nu, int M, int N){ +#define Fmu(A) Gimpl::CovShiftForward(Umu, mu, A) +#define Bmu(A) Gimpl::CovShiftBackward(Umu, mu, A) +#define Fnu(A) Gimpl::CovShiftForward(Unu, nu, A) +#define Bnu(A) Gimpl::CovShiftBackward(Unu, nu, A) +#define FmuI Gimpl::CovShiftIdentityForward(Umu, mu) +#define BmuI Gimpl::CovShiftIdentityBackward(Umu, mu) +#define FnuI Gimpl::CovShiftIdentityForward(Unu, nu) +#define BnuI Gimpl::CovShiftIdentityBackward(Unu, nu) + + //Upper right loop + GaugeMat tmp = BmuI; + for(int i=1;i(U, mu); + GaugeMat Unu = PeekIndex(U, nu); + if(M == N){ + GaugeMat F(Umu.Grid()); + CloverleafMxN(F, Umu, Unu, mu, nu, M, N); + FS = 0.125 * ( F - adj(F) ); + }else{ + //Average over both orientations + GaugeMat horizontal(Umu.Grid()), vertical(Umu.Grid()); + CloverleafMxN(horizontal, Umu, Unu, mu, nu, M, N); + CloverleafMxN(vertical, Umu, Unu, mu, nu, N, M); + FS = 0.0625 * ( horizontal - adj(horizontal) + vertical - adj(vertical) ); + } + } + + //Topological charge contribution from MxN Wilson loops + //cf https://arxiv.org/pdf/hep-lat/9701012.pdf Eq 6 + //output is the charge by timeslice: sum over timeslices to obtain the total + static std::vector TimesliceTopologicalChargeMxN(const GaugeLorentz &U, int M, int N){ + assert(Nd == 4); + std::vector > F(Nd,std::vector(Nd,nullptr)); + //Note F_numu = - F_munu + //hence we only need to loop over mu,nu,rho,sigma that aren't related by permuting mu,nu or rho,sigma + //Use nu > mu + for(int mu=0;mu Tq; + sliceSum(fsum, Tq, Nd-1); + + std::vector out(Tq.size()); + for(int t=0;t Tq = TimesliceTopologicalChargeMxN(U,M,N); + Real out(0); + for(int t=0;t > TimesliceTopologicalCharge5LiContributions(const GaugeLorentz &U){ + static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} }; + std::vector > out(5); + for(int i=0;i<5;i++){ + out[i] = TimesliceTopologicalChargeMxN(U,exts[i][0],exts[i][1]); + } + return out; + } + + static std::vector TopologicalCharge5LiContributions(const GaugeLorentz &U){ + static const int exts[5][2] = { {1,1}, {2,2}, {1,2}, {1,3}, {3,3} }; + std::vector out(5); + std::cout << GridLogMessage << "Computing topological charge" << std::endl; + for(int i=0;i<5;i++){ + out[i] = TopologicalChargeMxN(U,exts[i][0],exts[i][1]); + std::cout << GridLogMessage << exts[i][0] << "x" << exts[i][1] << " Wilson loop contribution " << out[i] << std::endl; + } + return out; + } + + //Compute the 5Li topological charge + static std::vector TimesliceTopologicalCharge5Li(const GaugeLorentz &U){ + std::vector > loops = TimesliceTopologicalCharge5LiContributions(U); + + double c5=1./20.; + double c4=1./5.-2.*c5; + double c3=(-64.+640.*c5)/45.; + double c2=(1-64.*c5)/9.; + double c1=(19.-55.*c5)/9.; + + int Lt = loops[0].size(); + std::vector out(Lt,0.); + for(int t=0;t Qt = TimesliceTopologicalCharge5Li(U); + Real Q = 0.; + for(int t=0;t::RandomGaugeTransform(RNG,U_GT,g); // Unit gauge + SU::RandomGaugeTransform(RNG,U_GT,g); // Unit gauge Field in_GT(&Grid); Field out_GT(&Grid); diff --git a/tests/core/Test_fft_gfix.cc b/tests/core/Test_fft_gfix.cc index 87dbc242..6d617e25 100644 --- a/tests/core/Test_fft_gfix.cc +++ b/tests/core/Test_fft_gfix.cc @@ -29,14 +29,10 @@ Author: Peter Boyle #include using namespace Grid; - ; -int main (int argc, char ** argv) -{ +template +void run(double alpha, bool do_fft_gfix){ std::vector seeds({1,2,3,4}); - - Grid_init(&argc,&argv); - int threads = GridThread::GetThreads(); Coordinate latt_size = GridDefaultLatt(); @@ -55,10 +51,7 @@ int main (int argc, char ** argv) FFT theFFT(&GRID); std::cout<::ColdConfiguration(pRNG,Umu); // Unit gauge Uorg=Umu; + + Real init_plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Initial plaquette "<< init_plaq << std::endl; + + //Apply a random gauge transformation to the unit gauge config Urnd=Umu; + SU::RandomGaugeTransform(pRNG,Urnd,g); - SU::RandomGaugeTransform(pRNG,Urnd,g); // Unit gauge - - Real plaq=WilsonLoops::avgPlaquette(Umu); - std::cout << " Initial plaquette "<::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false); + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform1,alpha,10000,1.0e-12, 1.0e-12,false); // Check the gauge xform matrices Utmp=Urnd; - SU::GaugeTransform(Utmp,xform1); + SU::GaugeTransform(Utmp,xform1); Utmp = Utmp - Umu; - std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl; + std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl; - plaq=WilsonLoops::avgPlaquette(Umu); - std::cout << " Final plaquette "<::avgPlaquette(Umu); + std::cout << " Final plaquette "<::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true); + + Utmp=Urnd; + SU::GaugeTransform(Utmp,xform2); + Utmp = Utmp - Umu; + std::cout << " Check the output gauge transformation matrices applied to the original field produce the xformed field "<< norm2(Utmp) << " (expect 0)" << std::endl; - std::cout<< "*****************************************************************" <::SteepestDescentGaugeFix(Umu,xform2,alpha,10000,1.0e-12, 1.0e-12,true); + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "<::GaugeTransform(Utmp,xform2); - Utmp = Utmp - Umu; - std::cout << " Norm Difference of xformed gauge "<< norm2(Utmp) << std::endl; + std::cout<< "******************************************************************************************" <::HotConfiguration(pRNG,Umu); - plaq=WilsonLoops::avgPlaquette(Umu); - std::cout << " Final plaquette "<::avgPlaquette(Umu); + std::cout << " Initial plaquette "<< init_plaq << std::endl; - std::cout<< "*****************************************************************" <::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,false); - SU::HotConfiguration(pRNG,Umu); // Unit gauge + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "<::avgPlaquette(Umu); - std::cout << " Initial plaquette "<::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true); + SU::HotConfiguration(pRNG,Umu); - plaq=WilsonLoops::avgPlaquette(Umu); - std::cout << " Final plaquette "<::avgPlaquette(Umu); + std::cout << " Initial plaquette "<< init_plaq << std::endl; - std::cout<< "*****************************************************************" <::SteepestDescentGaugeFix(Umu,alpha,10000,1.0e-12, 1.0e-12,true); + + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "<::HotConfiguration(pRNG,Umu); // Unit gauge + SU::HotConfiguration(pRNG,Umu); - plaq=WilsonLoops::avgPlaquette(Umu); - std::cout << " Initial plaquette "<::avgPlaquette(Umu); + std::cout << " Initial plaquette "<< init_plaq << std::endl; - FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir); + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,false,coulomb_dir); - std::cout << Umu<::avgPlaquette(Umu); + std::cout << " Final plaquette "<::avgPlaquette(Umu); - std::cout << " Final plaquette "<::HotConfiguration(pRNG,Umu); + + init_plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Initial plaquette "<< init_plaq << std::endl; + + FourierAcceleratedGaugeFixer::SteepestDescentGaugeFix(Umu,xform3,alpha,10000,1.0e-12, 1.0e-12,true,coulomb_dir); + + plaq=WilsonLoops::avgPlaquette(Umu); + std::cout << " Final plaquette "<> alpha; + } + } + + + if(gimpl == "periodic"){ + std::cout << GridLogMessage << "Using periodic boundary condition" << std::endl; + run(alpha, do_fft_gfix); + }else{ + std::vector conjdirs = {1,1,0,0}; //test with 2 conjugate dirs and 2 not + std::cout << GridLogMessage << "Using complex conjugate boundary conditions in dimensions "; + for(int i=0;i(alpha, do_fft_gfix); + } + Grid_finalize(); } diff --git a/tests/core/Test_gamma.cc b/tests/core/Test_gamma.cc index e52049fe..05f8c505 100644 --- a/tests/core/Test_gamma.cc +++ b/tests/core/Test_gamma.cc @@ -228,6 +228,59 @@ void checkGammaL(const Gamma::Algebra a, GridSerialRNG &rng) std::cout << std::endl; } +void checkChargeConjMatrix(){ + //Check the properties of the charge conjugation matrix + //In the Grid basis C = -\gamma^2 \gamma^4 + SpinMatrix C = testAlgebra[Gamma::Algebra::MinusGammaY] * testAlgebra[Gamma::Algebra::GammaT]; + SpinMatrix mC = -C; + SpinMatrix one = testAlgebra[Gamma::Algebra::Identity]; + + std::cout << "Testing properties of charge conjugation matrix C = -\\gamma^2 \\gamma^4 (in Grid's basis)" << std::endl; + + //C^T = -C + SpinMatrix Ct = transpose(C); + std::cout << GridLogMessage << "C^T=-C "; + test(Ct, mC); + std::cout << std::endl; + + //C^\dagger = -C + SpinMatrix Cdag = adj(C); + std::cout << GridLogMessage << "C^dag=-C "; + test(Cdag, mC); + std::cout << std::endl; + + //C^* = C + SpinMatrix Cstar = conjugate(C); + std::cout << GridLogMessage << "C^*=C "; + test(Cstar, C); + std::cout << std::endl; + + //C^{-1} = -C + SpinMatrix CinvC = mC * C; + std::cout << GridLogMessage << "C^{-1}=-C "; + test(CinvC, one); + std::cout << std::endl; + + // C^{-1} \gamma^\mu C = -[\gamma^\mu]^T + Gamma::Algebra gmu_a[4] = { Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ, Gamma::Algebra::GammaT }; + for(int mu=0;mu<4;mu++){ + SpinMatrix gmu = testAlgebra[gmu_a[mu]]; + SpinMatrix Cinv_gmu_C = mC * gmu * C; + SpinMatrix mgmu_T = -transpose(gmu); + std::cout << GridLogMessage << "C^{-1} \\gamma^" << mu << " C = -[\\gamma^" << mu << "]^T "; + test(Cinv_gmu_C, mgmu_T); + std::cout << std::endl; + } + + //[C, \gamma^5] = 0 + SpinMatrix Cg5 = C * testAlgebra[Gamma::Algebra::Gamma5]; + SpinMatrix g5C = testAlgebra[Gamma::Algebra::Gamma5] * C; + std::cout << GridLogMessage << "C \\gamma^5 = \\gamma^5 C"; + test(Cg5, g5C); + std::cout << std::endl; +} + + int main(int argc, char *argv[]) { Grid_init(&argc,&argv); @@ -270,6 +323,13 @@ int main(int argc, char *argv[]) { checkGammaL(i, sRNG); } + + std::cout << GridLogMessage << "======== Charge conjugation matrix check" << std::endl; + checkChargeConjMatrix(); + std::cout << GridLogMessage << std::endl; + + + Grid_finalize();