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; +}