diff --git a/lib/lattice/Lattice_base.h b/lib/lattice/Lattice_base.h index d97b1204..c91b7989 100644 --- a/lib/lattice/Lattice_base.h +++ b/lib/lattice/Lattice_base.h @@ -1,32 +1,33 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./lib/lattice/Lattice_base.h +Source file: ./lib/lattice/Lattice_base.h - Copyright (C) 2015 +Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle Author: paboyle - 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 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. +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. +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 */ +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ #ifndef GRID_LATTICE_BASE_H #define GRID_LATTICE_BASE_H @@ -255,6 +256,18 @@ PARALLEL_FOR_LOOP checkerboard=0; } + Lattice(const Lattice& r){ // copy constructor + _grid = r._grid; + checkerboard = r.checkerboard; + _odata.resize(_grid->oSites());// essential + PARALLEL_FOR_LOOP + for(int ss=0;ss<_grid->oSites();ss++){ + _odata[ss]=r._odata[ss]; + } + } + + + virtual ~Lattice(void) = default; template strong_inline Lattice & operator = (const sobj & r){ @@ -267,7 +280,7 @@ PARALLEL_FOR_LOOP template strong_inline Lattice & operator = (const Lattice & r){ this->checkerboard = r.checkerboard; conformable(*this,r); - std::cout< - class SmearedConfiguration { - public: - INHERIT_GIMPL_TYPES(Gimpl) ; +/*! + @brief Smeared configuration container - private: - const unsigned int smearingLevels; - Smear_Stout StoutSmearing; - std::vector SmearedSet; + It will behave like a configuration from the point of view of + the HMC update and integrators. + An "advanced configuration" object that can provide not only the + data to store the gauge configuration but also operations to manipulate + it, like smearing. - // Member functions - //==================================================================== - void fill_smearedSet(GaugeField& U){ - ThinLinks = &U; //attach the smearing routine to the field U + It stores a list of smeared configurations. +*/ +template +class SmearedConfiguration { + public: + INHERIT_GIMPL_TYPES(Gimpl); - //check the pointer is not null - if (ThinLinks==NULL) - std::cout << GridLogError << "[SmearedConfiguration] Error in ThinLinks pointer\n"; - - if (smearingLevels > 0){ - std::cout<< GridLogDebug << "[SmearedConfiguration] Filling SmearedSet\n"; - GaugeField previous_u(ThinLinks->_grid); + private: + const unsigned int smearingLevels; + Smear_Stout StoutSmearing; + std::vector SmearedSet; - previous_u = *ThinLinks; - for(int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl){ - StoutSmearing.smear(SmearedSet[smearLvl],previous_u); - previous_u = SmearedSet[smearLvl]; + // Member functions + //==================================================================== + void fill_smearedSet(GaugeField& U) { + ThinLinks = &U; // attach the smearing routine to the field U - // For debug purposes - RealD impl_plaq = WilsonLoops::avgPlaquette(previous_u); - std::cout<< GridLogDebug << "[SmearedConfiguration] Plaq: " << impl_plaq<< std::endl; + // check the pointer is not null + if (ThinLinks == NULL) + std::cout << GridLogError + << "[SmearedConfiguration] Error in ThinLinks pointer\n"; - } + if (smearingLevels > 0) { + std::cout << GridLogDebug + << "[SmearedConfiguration] Filling SmearedSet\n"; + GaugeField previous_u(ThinLinks->_grid); - } -} -//==================================================================== -GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, - const GaugeField& GaugeK) const{ - GridBase *grid = GaugeK._grid; - GaugeField C(grid), SigmaK(grid), iLambda(grid); - GaugeLinkField iLambda_mu(grid); - GaugeLinkField iQ(grid), e_iQ(grid); - GaugeLinkField SigmaKPrime_mu(grid); - GaugeLinkField GaugeKmu(grid), Cmu(grid); + previous_u = *ThinLinks; + for (int smearLvl = 0; smearLvl < smearingLevels; ++smearLvl) { + StoutSmearing.smear(SmearedSet[smearLvl], previous_u); + previous_u = SmearedSet[smearLvl]; - StoutSmearing.BaseSmear(C, GaugeK); - SigmaK = zero; - iLambda = zero; - - for (int mu = 0; mu < Nd; mu++){ - Cmu = peekLorentz( C,mu); - GaugeKmu = peekLorentz(GaugeK,mu); - SigmaKPrime_mu = peekLorentz(SigmaKPrime,mu); - iQ = Ta(Cmu*adj(GaugeKmu)); - set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu); - pokeLorentz(SigmaK, SigmaKPrime_mu*e_iQ + adj(Cmu)*iLambda_mu, mu); - pokeLorentz(iLambda, iLambda_mu, mu); - } - StoutSmearing.derivative(SigmaK, iLambda, GaugeK);// derivative of SmearBase - return SigmaK; -} - - - -/*! @brief Returns smeared configuration at level 'Level' */ -const GaugeField& get_smeared_conf(int Level) const{ - return SmearedSet[Level]; -} - - -//==================================================================== -void set_iLambda(GaugeLinkField& iLambda, - GaugeLinkField& e_iQ, - const GaugeLinkField& iQ, - const GaugeLinkField& Sigmap, - const GaugeLinkField& GaugeK)const{ - GridBase *grid = iQ._grid; - GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid); - GaugeLinkField unity(grid); - unity=1.0; - - LatticeComplex u(grid), w(grid); - LatticeComplex f0(grid), f1(grid), f2(grid); - LatticeComplex xi0(grid), xi1(grid), tmp(grid); - LatticeComplex u2(grid), w2(grid), cosw(grid); - LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid); - LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid); - LatticeComplex r22(grid), tr1(grid), tr2(grid); - LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), b22(grid); - LatticeComplex LatticeUnitComplex(grid); - - LatticeUnitComplex = 1.0; - - // Exponential - iQ2 = iQ * iQ; - iQ3 = iQ * iQ2; - StoutSmearing.set_uw(u,w,iQ2,iQ3); - StoutSmearing.set_fj(f0,f1,f2,u,w); - e_iQ = f0*unity + timesMinusI(f1) * iQ - f2 * iQ2; - - // Getting B1, B2, Gamma and Lambda - // simplify this part, reduntant calculations in set_fj - xi0 = StoutSmearing.func_xi0(w); - xi1 = StoutSmearing.func_xi1(w); - u2 = u * u; - w2 = w * w; - cosw = cos(w); - - emiu = cos(u) - timesI(sin(u)); - e2iu = cos(2.0*u) + timesI(sin(2.0*u)); - - r01 = (2.0*u + timesI(2.0*(u2-w2))) * e2iu - + emiu * ((16.0*u*cosw + 2.0*u*(3.0*u2+w2)*xi0) + - timesI(-8.0*u2*cosw + 2.0*(9.0*u2+w2)*xi0)); - - r11 = (2.0*LatticeUnitComplex + timesI(4.0*u))* e2iu - + emiu * ((-2.0*cosw + (3.0*u2-w2)*xi0) + - timesI((2.0*u*cosw + 6.0*u*xi0))); - - r21 = 2.0*timesI(e2iu) - + emiu * (-3.0*u*xi0 + timesI(cosw - 3.0*xi0)); - - - r02 = -2.0 * e2iu + emiu * (-8.0*u2*xi0 + - timesI(2.0*u*(cosw + xi0 + 3.0*u2*xi1))); - - r12 = emiu * (2.0*u*xi0 + timesI(-cosw - xi0 + 3.0*u2*xi1)); - - r22 = emiu * (xi0 - timesI(3.0*u*xi1)); - - fden = LatticeUnitComplex/(2.0*(9.0*u2-w2)*(9.0*u2-w2)); - - b10 = 2.0 * u * r01 + (3.0* u2 - w2)*r02 - (30.0 * u2 + 2.0 * w2)*f0; - b11 = 2.0 * u * r11 + (3.0* u2 - w2)*r12 - (30.0 * u2 + 2.0 * w2)*f1; - b12 = 2.0 * u * r21 + (3.0* u2 - w2)*r22 - (30.0 * u2 + 2.0 * w2)*f2; - - b20 = r01 - (3.0*u)*r02 - (24.0*u)*f0; - b21 = r11 - (3.0*u)*r12 - (24.0*u)*f1; - b22 = r21 - (3.0*u)*r22 - (24.0*u)*f2; - - b10 *= fden; - b11 *= fden; - b12 *= fden; - b20 *= fden; - b21 *= fden; - b22 *= fden; - - - B1 = b10*unity + timesMinusI(b11) * iQ - b12 * iQ2; - B2 = b20*unity + timesMinusI(b21) * iQ - b22 * iQ2; - USigmap = GaugeK * Sigmap; - - tr1 = trace(USigmap*B1); - tr2 = trace(USigmap*B2); - - GaugeLinkField QUS = iQ * USigmap; - GaugeLinkField USQ = USigmap * iQ; - - GaugeLinkField iGamma = tr1 * iQ - timesI(tr2) * iQ2 + - timesI(f1) * USigmap + f2 * QUS + f2 * USQ; - - iLambda = Ta(iGamma); - -} - -//==================================================================== -public: - GaugeField* ThinLinks; /*!< @brief Pointer to the thin - links configuration */ - - /*! @brief Standard constructor */ - SmearedConfiguration(GridCartesian * UGrid, - unsigned int Nsmear, - Smear_Stout& Stout): - smearingLevels(Nsmear), - StoutSmearing(Stout), - ThinLinks(NULL){ - for (unsigned int i=0; i< smearingLevels; ++i) - SmearedSet.push_back(*(new GaugeField(UGrid))); + // For debug purposes + RealD impl_plaq = WilsonLoops::avgPlaquette(previous_u); + std::cout << GridLogDebug + << "[SmearedConfiguration] Plaq: " << impl_plaq << std::endl; + } } + } + //==================================================================== + GaugeField AnalyticSmearedForce(const GaugeField& SigmaKPrime, + const GaugeField& GaugeK) const { + GridBase* grid = GaugeK._grid; + GaugeField C(grid), SigmaK(grid), iLambda(grid); + GaugeLinkField iLambda_mu(grid); + GaugeLinkField iQ(grid), e_iQ(grid); + GaugeLinkField SigmaKPrime_mu(grid); + GaugeLinkField GaugeKmu(grid), Cmu(grid); - /*! For just thin links */ - SmearedConfiguration(): - smearingLevels(0), - StoutSmearing(), - SmearedSet(), - ThinLinks(NULL){} + StoutSmearing.BaseSmear(C, GaugeK); + SigmaK = zero; + iLambda = zero; + for (int mu = 0; mu < Nd; mu++) { + Cmu = peekLorentz(C, mu); + GaugeKmu = peekLorentz(GaugeK, mu); + SigmaKPrime_mu = peekLorentz(SigmaKPrime, mu); + iQ = Ta(Cmu * adj(GaugeKmu)); + set_iLambda(iLambda_mu, e_iQ, iQ, SigmaKPrime_mu, GaugeKmu); + pokeLorentz(SigmaK, SigmaKPrime_mu * e_iQ + adj(Cmu) * iLambda_mu, mu); + pokeLorentz(iLambda, iLambda_mu, mu); + } + StoutSmearing.derivative(SigmaK, iLambda, + GaugeK); // derivative of SmearBase + return SigmaK; + } - // attach the smeared routines to the thin links U and fill the smeared set - void set_GaugeField(GaugeField& U){ fill_smearedSet(U);} + /*! @brief Returns smeared configuration at level 'Level' */ + const GaugeField& get_smeared_conf(int Level) const { + return SmearedSet[Level]; + } -//==================================================================== - void smeared_force(GaugeField& SigmaTilde) const{ + //==================================================================== + void set_iLambda(GaugeLinkField& iLambda, GaugeLinkField& e_iQ, + const GaugeLinkField& iQ, const GaugeLinkField& Sigmap, + const GaugeLinkField& GaugeK) const { + GridBase* grid = iQ._grid; + GaugeLinkField iQ2(grid), iQ3(grid), B1(grid), B2(grid), USigmap(grid); + GaugeLinkField unity(grid); + unity = 1.0; - if (smearingLevels > 0){ - GaugeField force(SigmaTilde._grid); - GaugeLinkField tmp_mu(SigmaTilde._grid); - force = SigmaTilde;//actually = U*SigmaTilde + LatticeComplex u(grid), w(grid); + LatticeComplex f0(grid), f1(grid), f2(grid); + LatticeComplex xi0(grid), xi1(grid), tmp(grid); + LatticeComplex u2(grid), w2(grid), cosw(grid); + LatticeComplex emiu(grid), e2iu(grid), qt(grid), fden(grid); + LatticeComplex r01(grid), r11(grid), r21(grid), r02(grid), r12(grid); + LatticeComplex r22(grid), tr1(grid), tr2(grid); + LatticeComplex b10(grid), b11(grid), b12(grid), b20(grid), b21(grid), + b22(grid); + LatticeComplex LatticeUnitComplex(grid); - for (int mu = 0; mu < Nd; mu++){ - // to get just SigmaTilde - tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels-1], mu)) * peekLorentz(force,mu); - pokeLorentz(force, tmp_mu, mu); - } + LatticeUnitComplex = 1.0; - for(int ismr = smearingLevels - 1; ismr > 0; --ismr) - force = AnalyticSmearedForce(force,get_smeared_conf(ismr-1)); + // Exponential + iQ2 = iQ * iQ; + iQ3 = iQ * iQ2; + StoutSmearing.set_uw(u, w, iQ2, iQ3); + StoutSmearing.set_fj(f0, f1, f2, u, w); + e_iQ = f0 * unity + timesMinusI(f1) * iQ - f2 * iQ2; - force = AnalyticSmearedForce(force,*ThinLinks); - - for (int mu = 0; mu < Nd; mu++){ - tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); - pokeLorentz(SigmaTilde, tmp_mu, mu); - } - }// if smearingLevels = 0 do nothing -} - //==================================================================== + // Getting B1, B2, Gamma and Lambda + // simplify this part, reduntant calculations in set_fj + xi0 = StoutSmearing.func_xi0(w); + xi1 = StoutSmearing.func_xi1(w); + u2 = u * u; + w2 = w * w; + cosw = cos(w); + emiu = cos(u) - timesI(sin(u)); + e2iu = cos(2.0 * u) + timesI(sin(2.0 * u)); -GaugeField& get_SmearedU(){ - return SmearedSet[smearingLevels-1]; -} + r01 = (2.0 * u + timesI(2.0 * (u2 - w2))) * e2iu + + emiu * ((16.0 * u * cosw + 2.0 * u * (3.0 * u2 + w2) * xi0) + + timesI(-8.0 * u2 * cosw + 2.0 * (9.0 * u2 + w2) * xi0)); -GaugeField& get_U(bool smeared=false) { - // get the config, thin links by default - if (smeared){ - if (smearingLevels){ - RealD impl_plaq = WilsonLoops::avgPlaquette(SmearedSet[smearingLevels-1]); - std::cout<< GridLogDebug << "getting Usmr Plaq: " << impl_plaq<< std::endl; - return get_SmearedU(); + r11 = (2.0 * LatticeUnitComplex + timesI(4.0 * u)) * e2iu + + emiu * ((-2.0 * cosw + (3.0 * u2 - w2) * xi0) + + timesI((2.0 * u * cosw + 6.0 * u * xi0))); - } - else { - RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); - std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; - return *ThinLinks; - } - } - else{ - RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); - std::cout<< GridLogDebug << "getting Thin Plaq: " << impl_plaq<< std::endl; - return *ThinLinks;} - } + r21 = + 2.0 * timesI(e2iu) + emiu * (-3.0 * u * xi0 + timesI(cosw - 3.0 * xi0)); + r02 = -2.0 * e2iu + + emiu * (-8.0 * u2 * xi0 + + timesI(2.0 * u * (cosw + xi0 + 3.0 * u2 * xi1))); + + r12 = emiu * (2.0 * u * xi0 + timesI(-cosw - xi0 + 3.0 * u2 * xi1)); + + r22 = emiu * (xi0 - timesI(3.0 * u * xi1)); + + fden = LatticeUnitComplex / (2.0 * (9.0 * u2 - w2) * (9.0 * u2 - w2)); + + b10 = 2.0 * u * r01 + (3.0 * u2 - w2) * r02 - (30.0 * u2 + 2.0 * w2) * f0; + b11 = 2.0 * u * r11 + (3.0 * u2 - w2) * r12 - (30.0 * u2 + 2.0 * w2) * f1; + b12 = 2.0 * u * r21 + (3.0 * u2 - w2) * r22 - (30.0 * u2 + 2.0 * w2) * f2; + + b20 = r01 - (3.0 * u) * r02 - (24.0 * u) * f0; + b21 = r11 - (3.0 * u) * r12 - (24.0 * u) * f1; + b22 = r21 - (3.0 * u) * r22 - (24.0 * u) * f2; + + b10 *= fden; + b11 *= fden; + b12 *= fden; + b20 *= fden; + b21 *= fden; + b22 *= fden; + + B1 = b10 * unity + timesMinusI(b11) * iQ - b12 * iQ2; + B2 = b20 * unity + timesMinusI(b21) * iQ - b22 * iQ2; + USigmap = GaugeK * Sigmap; + + tr1 = trace(USigmap * B1); + tr2 = trace(USigmap * B2); + + GaugeLinkField QUS = iQ * USigmap; + GaugeLinkField USQ = USigmap * iQ; + + GaugeLinkField iGamma = tr1 * iQ - timesI(tr2) * iQ2 + + timesI(f1) * USigmap + f2 * QUS + f2 * USQ; + + iLambda = Ta(iGamma); + } + + //==================================================================== + public: + GaugeField* + ThinLinks; /*!< @brief Pointer to the thin + links configuration */ + + /*! @brief Standard constructor */ + SmearedConfiguration(GridCartesian* UGrid, unsigned int Nsmear, + Smear_Stout& Stout) + : smearingLevels(Nsmear), StoutSmearing(Stout), ThinLinks(NULL) { + for (unsigned int i = 0; i < smearingLevels; ++i) + SmearedSet.push_back(*(new GaugeField(UGrid))); + } + + /*! For just thin links */ + SmearedConfiguration() + : smearingLevels(0), StoutSmearing(), SmearedSet(), ThinLinks(NULL) {} + + // attach the smeared routines to the thin links U and fill the smeared set + void set_GaugeField(GaugeField& U) { fill_smearedSet(U); } + + //==================================================================== + void smeared_force(GaugeField& SigmaTilde) const { + if (smearingLevels > 0) { + GaugeField force = SigmaTilde; // actually = U*SigmaTilde + GaugeLinkField tmp_mu(SigmaTilde._grid); + + for (int mu = 0; mu < Nd; mu++) { + // to get just SigmaTilde + tmp_mu = adj(peekLorentz(SmearedSet[smearingLevels - 1], mu)) * + peekLorentz(force, mu); + pokeLorentz(force, tmp_mu, mu); + } + + for (int ismr = smearingLevels - 1; ismr > 0; --ismr) + force = AnalyticSmearedForce(force, get_smeared_conf(ismr - 1)); + + force = AnalyticSmearedForce(force, *ThinLinks); + + for (int mu = 0; mu < Nd; mu++) { + tmp_mu = peekLorentz(*ThinLinks, mu) * peekLorentz(force, mu); + pokeLorentz(SigmaTilde, tmp_mu, mu); + } + } // if smearingLevels = 0 do nothing + } + //==================================================================== + + GaugeField& get_SmearedU() { return SmearedSet[smearingLevels - 1]; } + + GaugeField& get_U(bool smeared = false) { + // get the config, thin links by default + if (smeared) { + if (smearingLevels) { + RealD impl_plaq = + WilsonLoops::avgPlaquette(SmearedSet[smearingLevels - 1]); + std::cout << GridLogDebug << "getting Usmr Plaq: " << impl_plaq + << std::endl; + return get_SmearedU(); + + } else { + RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); + std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq + << std::endl; + return *ThinLinks; + } + } else { + RealD impl_plaq = WilsonLoops::avgPlaquette(*ThinLinks); + std::cout << GridLogDebug << "getting Thin Plaq: " << impl_plaq + << std::endl; + return *ThinLinks; + } + } }; - - } - } - - - - - #endif diff --git a/scripts/copyright b/scripts/copyright index 76ed6603..92772f16 100755 --- a/scripts/copyright +++ b/scripts/copyright @@ -5,13 +5,13 @@ while (( "$#" )); do echo $1 cat > message <> message cat >> message < tmp.fil diff --git a/tests/Test_main.cc b/tests/Test_main.cc index 3cbe5229..7eb0107b 100644 --- a/tests/Test_main.cc +++ b/tests/Test_main.cc @@ -1,620 +1,666 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/Test_main.cc +Source file: ./tests/Test_main.cc - Copyright (C) 2015 +Copyright (C) 2015 Author: Azusa Yamaguchi Author: Peter Boyle Author: neo Author: paboyle - 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 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. +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. +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 */ +See the full license in the file "LICENSE" in the top level distribution +directory +*************************************************************************************/ +/* END LEGAL */ #include "Grid.h" - using namespace std; using namespace Grid; using namespace Grid::QCD; /* - Grid_main.cc(232): error: no suitable user-defined conversion from "Grid::iScalar, 4>>" to "const Grid::iScalar>>" exists + Grid_main.cc(232): error: no suitable user-defined conversion from +"Grid::iScalar, 4>>" to "const +Grid::iScalar>>" exists c_m = peekIdiot(scm,1,2); */ -template auto peekIdiot(const vobj &rhs,int i,int j) -> decltype(peekIndex<2>(rhs,0,0)) -{ - return peekIndex<2>(rhs,i,j); +template +auto peekIdiot(const vobj &rhs, int i, int j) + -> decltype(peekIndex<2>(rhs, 0, 0)) { + return peekIndex<2>(rhs, i, j); } -template auto peekDumKopf(const vobj &rhs,int i,int j) -> decltype(peekIndex<3>(rhs,0,0)) -{ - return peekIndex<3>(rhs,i,j); +template +auto peekDumKopf(const vobj &rhs, int i, int j) + -> decltype(peekIndex<3>(rhs, 0, 0)) { + return peekIndex<3>(rhs, i, j); } -template auto peekDumKopf(const vobj &rhs,int i) -> decltype(peekIndex<3>(rhs,0)) -{ - return peekIndex<3>(rhs,i); +template +auto peekDumKopf(const vobj &rhs, int i) -> decltype(peekIndex<3>(rhs, 0)) { + return peekIndex<3>(rhs, i); } -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); +int main(int argc, char **argv) { + Grid_init(&argc, &argv); - std::vector latt_size = GridDefaultLatt(); - std::vector simd_layout = GridDefaultSimd(4,vComplex::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); + std::vector latt_size = GridDefaultLatt(); + std::vector simd_layout = GridDefaultSimd(4, vComplex::Nsimd()); + std::vector mpi_layout = GridDefaultMpi(); latt_size.resize(4); #ifdef AVX512 - for(int omp=128;omp<236;omp+=16){ + for (int omp = 128; omp < 236; omp += 16) { #else - for(int omp=1;omp<2;omp*=20){ + for (int omp = 1; omp < 2; omp *= 20) { #endif #ifdef OMP - omp_set_num_threads(omp); -#endif - - for(int lat=8;lat<=16;lat+=40){ - - std::cout << "Lat "< saved; - SerialRNG.GetState(saved,0); - SerialRNG1.SetState(saved,0); - - RealD dd1,dd2; - - std::cout << "Testing RNG state save restore"< iGammaFive; - ColourMatrix cmat; - - random(FineRNG,Foo); - gaussian(FineRNG,Bar); - random(FineRNG,scFoo); - random(FineRNG,scBar); - - random(FineRNG,cMat); - random(FineRNG,sMat); - random(FineRNG,scMat); - random(FineRNG,lcMat); - random(FineRNG,cVec); - random(FineRNG,sVec); - random(FineRNG,scVec); - - - fflush(stdout); - - TComplex tr = trace(cmat); - - - cVec = cMat * cVec; // LatticeColourVector = LatticeColourMatrix * LatticeColourVector - sVec = sMat * sVec; // LatticeSpinVector = LatticeSpinMatrix * LatticeSpinVector - scVec= scMat * scVec;// LatticeSpinColourVector = LatticeSpinColourMatrix * LatticeSpinColourVector - scVec= cMat * scVec; // LatticeSpinColourVector = LatticeColourMatrix * LatticeSpinColourVector - scVec= sMat * scVec; // LatticeSpinColourVector = LatticeSpinMatrix * LatticeSpinColourVector - - cMat = outerProduct(cVec,cVec); - scalar = localInnerProduct(cVec,cVec); - - cMat = Ta(cMat); //traceless antihermitian - - - scalar += scalar; - scalar -= scalar; - scalar *= scalar; - add(scalar,scalar,scalar); - sub(scalar,scalar,scalar); - mult(scalar,scalar,scalar); - - mac(scalar,scalar,scalar); - scalar = scalar+scalar; - scalar = scalar-scalar; - scalar = scalar*scalar; - - scalar=outerProduct(scalar,scalar); - - scalar=adj(scalar); - - // rscalar=real(scalar); - // iscalar=imag(scalar); - // scalar =cmplx(rscalar,iscalar); - PokeIndex(cVec,scalar,1); - - - scalar=transpose(scalar); - scalar=TransposeIndex(scalar); - scalar=TraceIndex(scalar); - scalar=PeekIndex(cVec,0); - - scalar=trace(scalar); - scalar=localInnerProduct(cVec,cVec); - scalar=localNorm2(cVec); - -// -=,+=,*=,() -// add,+,sub,-,mult,mac,* -// adj,conjugate -// real,imag -// transpose,transposeIndex -// trace,traceIndex -// peekIndex -// innerProduct,outerProduct, -// localNorm2 -// localInnerProduct - - - scMat = sMat*scMat; // LatticeSpinColourMatrix = LatticeSpinMatrix * LatticeSpinColourMatrix - - - /////////////////////// - // Non-lattice (const objects) * Lattice - ColourMatrix cm; - SpinColourMatrix scm; - vSpinColourMatrix vscm; - Complex cplx(1.0); - Integer myint=1; - double mydouble=1.0; - - // vSpinColourMatrix vscm; - scMat = cMat*scMat; - scm = cm * scm; // SpinColourMatrix = ColourMatrix * SpinColourMatrix - scm = scm *cm; // SpinColourMatrix = SpinColourMartix * ColourMatrix - scm = GammaFive * scm ; // SpinColourMatrix = SpinMatrix * SpinColourMatrix - scm = scm* GammaFive ; // SpinColourMatrix = SpinColourMatrix * SpinMatrix - - scm = scm*cplx; - vscm = vscm*cplx; - scMat = scMat*cplx; - - scm = cplx*scm; - vscm = cplx*vscm; - scMat = cplx*scMat; - scm = myint*scm; - vscm = myint*vscm; - scMat = scMat*myint; - - scm = scm*mydouble; - vscm = vscm*mydouble; - scMat = scMat*mydouble; - scMat = mydouble*scMat; - cMat = mydouble*cMat; - - sMat = adj(sMat); // LatticeSpinMatrix adjoint - sMat = iGammaFive*sMat; // SpinMatrix * LatticeSpinMatrix - sMat = GammaFive*sMat; // SpinMatrix * LatticeSpinMatrix - scMat= adj(scMat); - cMat= adj(cMat); - cm=adj(cm); - scm=adj(scm); - scm=transpose(scm); - scm=transposeIndex<1>(scm); - - - random(SerialRNG, cm); - std::cout< mysite {0,0,0,0}; - random(FineRNG,cMat); - cMat = Ta(cMat); - peekSite(cm, cMat, mysite); - std::cout<::traceIndex(sc_m); // Map to traceColour - c_m = TensorIndexRecursion::traceIndex(sc_m); // map to traceSpin - - c = TensorIndexRecursion::traceIndex(s_m); - c = TensorIndexRecursion::traceIndex(c_m); - - s_m = TensorIndexRecursion::peekIndex(scm,0,0); - c_m = TensorIndexRecursion::peekIndex(scm,1,2); - // c_m = peekSpin(scm,1,2); - // c_m = peekIdiot(scm,1,2); - - printf("c. Level %d\n",c_m.TensorLevel); - printf("c. Level %d\n",c_m().TensorLevel); - printf("c. Level %d\n",c_m()().TensorLevel); - - c_m()() = scm()(0,0); //ColourComponents of CM <= ColourComponents of SpinColourMatrix - scm()(1,1) = cm()(); //ColourComponents of CM <= ColourComponents of SpinColourMatrix - c = scm()(1,1)(1,2); - scm()(1,1)(2,1) = c; - - // pokeIndex (c_m,c,0,0); - } - - FooBar = Bar; - /* - { - std::vector coor(4); - for(int d=0;d<4;d++) coor[d] = 0; - peekSite(cmat,Foo,coor); - Foo = zero; - pokeSite(cmat,Foo,coor); - } - random(Foo); - */ - lex_sites(Foo); - - - Integer mm[4]; - mm[0]=1; - mm[1]=Fine._rdimensions[0]; - mm[2]=Fine._ldimensions[0]*Fine._ldimensions[1]; - mm[3]=Fine._ldimensions[0]*Fine._ldimensions[1]*Fine._ldimensions[2]; - - LatticeInteger lex(&Fine); - lex=zero; - for(int d=0;d<4;d++){ - LatticeInteger coor(&Fine); - LatticeCoordinate(coor,d); - lex = lex + coor*mm[d]; - - } - - - - // Bar = zero; - // Bar = where(lex coor(4); - for(coor[3]=0;coor[3] saved; + SerialRNG.GetState(saved, 0); + SerialRNG1.SetState(saved, 0); + + RealD dd1, dd2; + + std::cout << "Testing RNG state save restore" << std::endl; + for (int i = 0; i < 10; i++) { + random(SerialRNG, dd1); + random(SerialRNG1, dd2); + std::cout << "i " << i << " " << dd1 << " " << dd2 << std::endl; + } + LatticeColourMatrix Foo(&Fine); + LatticeColourMatrix Bar(&Fine); + + LatticeSpinColourMatrix scFoo(&Fine); + LatticeSpinColourMatrix scBar(&Fine); + + LatticeColourMatrix Shifted(&Fine); + LatticeColourMatrix ShiftedCheck(&Fine); + LatticeColourMatrix rShifted(&rbFine); + LatticeColourMatrix bShifted(&rbFine); + + LatticeColourMatrix rFoo(&rbFine); + LatticeColourMatrix bFoo(&rbFine); + + LatticeColourMatrix FooBar(&Fine); + LatticeSpinColourMatrix scFooBar(&Fine); + + LatticeColourVector cVec(&Fine); + LatticeSpinVector sVec(&Fine); + LatticeSpinColourVector scVec(&Fine); + + LatticeColourMatrix cMat(&Fine); + LatticeSpinMatrix sMat(&Fine); + LatticeSpinColourMatrix scMat(&Fine); + + LatticeLorentzColourMatrix lcMat(&Fine); + + LatticeComplex scalar(&Fine); + LatticeReal rscalar(&Fine); + LatticeReal iscalar(&Fine); + + SpinMatrix GammaFive; + iSpinMatrix iGammaFive; + ColourMatrix cmat; + + random(FineRNG, Foo); + gaussian(FineRNG, Bar); + random(FineRNG, scFoo); + random(FineRNG, scBar); + + random(FineRNG, cMat); + random(FineRNG, sMat); + random(FineRNG, scMat); + random(FineRNG, lcMat); + random(FineRNG, cVec); + random(FineRNG, sVec); + random(FineRNG, scVec); + + fflush(stdout); + + LatticeColourMatrix newFoo = Foo; + // confirm correctness of copy constructor + Bar = Foo - newFoo; + std::cout << "Copy constructor diff check: "; + double test_cc = norm2(Bar); + if (test_cc < 1e-5){ + std::cout << "OK\n"; + } + else{ + std::cout << "fail\n"; + abort(); + } + + TComplex tr = trace(cmat); + + cVec = cMat * cVec; // LatticeColourVector = LatticeColourMatrix + // * LatticeColourVector + sVec = sMat * sVec; // LatticeSpinVector = LatticeSpinMatrix + // * LatticeSpinVector + scVec = scMat * scVec; // LatticeSpinColourVector = + // LatticeSpinColourMatrix * + // LatticeSpinColourVector + scVec = cMat * scVec; // LatticeSpinColourVector = LatticeColourMatrix + // * LatticeSpinColourVector + scVec = sMat * scVec; // LatticeSpinColourVector = LatticeSpinMatrix + // * LatticeSpinColourVector + + cMat = outerProduct(cVec, cVec); + scalar = localInnerProduct(cVec, cVec); + + cMat = Ta(cMat); // traceless antihermitian + + scalar += scalar; + scalar -= scalar; + scalar *= scalar; + add(scalar, scalar, scalar); + sub(scalar, scalar, scalar); + mult(scalar, scalar, scalar); + + mac(scalar, scalar, scalar); + scalar = scalar + scalar; + scalar = scalar - scalar; + scalar = scalar * scalar; + + scalar = outerProduct(scalar, scalar); + + scalar = adj(scalar); + + // rscalar=real(scalar); + // iscalar=imag(scalar); + // scalar =cmplx(rscalar,iscalar); + PokeIndex(cVec, scalar, 1); + + scalar = transpose(scalar); + scalar = TransposeIndex(scalar); + scalar = TraceIndex(scalar); + scalar = PeekIndex(cVec, 0); + + scalar = trace(scalar); + scalar = localInnerProduct(cVec, cVec); + scalar = localNorm2(cVec); + + // -=,+=,*=,() + // add,+,sub,-,mult,mac,* + // adj,conjugate + // real,imag + // transpose,transposeIndex + // trace,traceIndex + // peekIndex + // innerProduct,outerProduct, + // localNorm2 + // localInnerProduct + + scMat = sMat * scMat; // LatticeSpinColourMatrix = LatticeSpinMatrix + // * LatticeSpinColourMatrix + + /////////////////////// + // Non-lattice (const objects) * Lattice + ColourMatrix cm; + SpinColourMatrix scm; + vSpinColourMatrix vscm; + Complex cplx(1.0); + Integer myint = 1; + double mydouble = 1.0; + + // vSpinColourMatrix vscm; + scMat = cMat * scMat; + scm = + cm * scm; // SpinColourMatrix = ColourMatrix * SpinColourMatrix + scm = scm * cm; // SpinColourMatrix = SpinColourMartix * ColourMatrix + scm = GammaFive * + scm; // SpinColourMatrix = SpinMatrix * SpinColourMatrix + scm = + scm * GammaFive; // SpinColourMatrix = SpinColourMatrix * SpinMatrix + + scm = scm * cplx; + vscm = vscm * cplx; + scMat = scMat * cplx; + + scm = cplx * scm; + vscm = cplx * vscm; + scMat = cplx * scMat; + scm = myint * scm; + vscm = myint * vscm; + scMat = scMat * myint; + + scm = scm * mydouble; + vscm = vscm * mydouble; + scMat = scMat * mydouble; + scMat = mydouble * scMat; + cMat = mydouble * cMat; + + sMat = adj(sMat); // LatticeSpinMatrix adjoint + sMat = iGammaFive * sMat; // SpinMatrix * LatticeSpinMatrix + sMat = GammaFive * sMat; // SpinMatrix * LatticeSpinMatrix + scMat = adj(scMat); + cMat = adj(cMat); + cm = adj(cm); + scm = adj(scm); + scm = transpose(scm); + scm = transposeIndex<1>(scm); + + random(SerialRNG, cm); + std::cout << GridLogMessage << cm << std::endl; + + cm = Ta(cm); + TComplex tracecm = trace(cm); + std::cout << GridLogMessage << cm << std::endl; + + cm = Exponentiate(cm, 2.0, 12); + std::cout << GridLogMessage << cm << " " << std::endl; + Complex det = Determinant(cm); + std::cout << GridLogMessage << "determinant: " << det << std::endl; + std::cout << GridLogMessage << "norm: " << norm2(cm) << std::endl; + + cm = ProjectOnGroup(cm); + std::cout << GridLogMessage << cm << " " << std::endl; + std::cout << GridLogMessage << "norm: " << norm2(cm) << std::endl; + cm = ProjectOnGroup(cm); + std::cout << GridLogMessage << cm << " " << std::endl; + std::cout << GridLogMessage << "norm: " << norm2(cm) << std::endl; + + // det = Determinant(cm); + // std::cout< mysite{0, 0, 0, 0}; + random(FineRNG, cMat); + cMat = Ta(cMat); + peekSite(cm, cMat, mysite); + std::cout << GridLogMessage << cm << " " << std::endl; + cm = Exponentiate(cm, 1.0, 12); + std::cout << GridLogMessage << cm << " " << std::endl; + std::cout << GridLogMessage << "norm: " << norm2(cm) << std::endl; + + std::cout << GridLogMessage << "norm cMmat : " << norm2(cMat) + << std::endl; + cMat = expMat(cMat, ComplexD(1.0, 0.0)); + std::cout << GridLogMessage << "norm expMat: " << norm2(cMat) + << std::endl; + peekSite(cm, cMat, mysite); + std::cout << GridLogMessage << cm << " " << std::endl; + std::cout << GridLogMessage << "determinant: " << Determinant(cm) + << std::endl; + std::cout << GridLogMessage << "norm: " << norm2(cm) << std::endl; + + // LatticeComplex trlcMat(&Fine); + // trlcMat = trace(lcMat); // Trace involving iVector - now generates + // error + + { // Peek-ology and Poke-ology, with a little app-ology + Complex c; + ColourMatrix c_m; + SpinMatrix s_m; + SpinColourMatrix sc_m; + + s_m = TensorIndexRecursion::traceIndex( + sc_m); // Map to traceColour + c_m = TensorIndexRecursion::traceIndex( + sc_m); // map to traceSpin + + c = TensorIndexRecursion::traceIndex(s_m); + c = TensorIndexRecursion::traceIndex(c_m); + + s_m = TensorIndexRecursion::peekIndex(scm, 0, 0); + c_m = TensorIndexRecursion::peekIndex(scm, 1, 2); + // c_m = peekSpin(scm,1,2); + // c_m = peekIdiot(scm,1,2); + + printf("c. Level %d\n", c_m.TensorLevel); + printf("c. Level %d\n", c_m().TensorLevel); + printf("c. Level %d\n", c_m()().TensorLevel); + + c_m()() = scm()(0, 0); // ColourComponents of CM <= ColourComponents of + // SpinColourMatrix + scm()(1, 1) = cm()(); // ColourComponents of CM <= ColourComponents of + // SpinColourMatrix + c = scm()(1, 1)(1, 2); + scm()(1, 1)(2, 1) = c; + + // pokeIndex (c_m,c,0,0); + } + + FooBar = Bar; + /* + { + std::vector coor(4); + for(int d=0;d<4;d++) coor[d] = 0; + peekSite(cmat,Foo,coor); + Foo = zero; + pokeSite(cmat,Foo,coor); + } + random(Foo); + */ + lex_sites(Foo); + + Integer mm[4]; + mm[0] = 1; + mm[1] = Fine._rdimensions[0]; + mm[2] = Fine._ldimensions[0] * Fine._ldimensions[1]; + mm[3] = + Fine._ldimensions[0] * Fine._ldimensions[1] * Fine._ldimensions[2]; + + LatticeInteger lex(&Fine); + lex = zero; + for (int d = 0; d < 4; d++) { + LatticeInteger coor(&Fine); + LatticeCoordinate(coor, d); + lex = lex + coor * mm[d]; + } + + // Bar = zero; + // Bar = where(lex coor(4); + for (coor[3] = 0; coor[3] < latt_size[3] / mpi_layout[3]; coor[3]++) { + for (coor[2] = 0; coor[2] < latt_size[2] / mpi_layout[2]; coor[2]++) { + for (coor[1] = 0; coor[1] < latt_size[1] / mpi_layout[1]; + coor[1]++) { + for (coor[0] = 0; coor[0] < latt_size[0] / mpi_layout[0]; + coor[0]++) { + ColourMatrix bar; + peekSite(bar, Bar, coor); + for (int r = 0; r < 3; r++) { + for (int c = 0; c < 3; c++) { + // cout<<"bar "<black - rShifted = Cshift(bFoo,dir,shift); // Shift black->red - - ShiftedCheck=zero; - setCheckerboard(ShiftedCheck,bShifted); // Put them all together - setCheckerboard(ShiftedCheck,rShifted); // and check the results (later) - - // Check results - std::vector coor(4); - for(coor[3]=0;coor[3] diff; - - std::vector shiftcoor = coor; - shiftcoor[dir]=(shiftcoor[dir]+shift+latt_size[dir])%(latt_size[dir]/mpi_layout[dir]); + if (Fine.IsBoss()) { + printf("Cshift Mult: NumThread %d , Lattice size %d , %f us per call\n", + omp, lat, (t1 - t0) / ncall); + printf("Cshift Mult: NumThread %d , Lattice size %d , %f Mflop/s\n", + omp, lat, flops / (t1 - t0)); + printf("Cshift Mult: NumThread %d , Lattice size %d , %f MB/s\n", omp, + lat, bytes / (t1 - t0)); + } + // pickCheckerboard(0,rFoo,FooBar); + // pickCheckerboard(1,bFoo,FooBar); + // setCheckerboard(FooBar,rFoo); + // setCheckerboard(FooBar,bFoo); - std::vector rl(4); - for(int dd=0;dd<4;dd++){ - rl[dd] = latt_size[dd]/simd_layout[dd]/mpi_layout[dd]; - } - int lex = coor[0]%rl[0] - + (coor[1]%rl[1])*rl[0] - + (coor[2]%rl[2])*rl[0]*rl[1] - + (coor[3]%rl[3])*rl[0]*rl[1]*rl[2]; - lex += - +1000*(coor[0]/rl[0]) - +1000*(coor[1]/rl[1])*simd_layout[0] - +1000*(coor[2]/rl[2])*simd_layout[0]*simd_layout[1] - +1000*(coor[3]/rl[3])*simd_layout[0]*simd_layout[1]*simd_layout[2]; + double nrm = 0; - int lex_coor = shiftcoor[0]%rl[0] - + (shiftcoor[1]%rl[1])*rl[0] - + (shiftcoor[2]%rl[2])*rl[0]*rl[1] - + (shiftcoor[3]%rl[3])*rl[0]*rl[1]*rl[2]; - lex_coor += - +1000*(shiftcoor[0]/rl[0]) - +1000*(shiftcoor[1]/rl[1])*simd_layout[0] - +1000*(shiftcoor[2]/rl[2])*simd_layout[0]*simd_layout[1] - +1000*(shiftcoor[3]/rl[3])*simd_layout[0]*simd_layout[1]*simd_layout[2]; + LatticeColourMatrix deriv(&Fine); + double half = 0.5; + deriv = 0.5 * Cshift(Foo, 0, 1) - 0.5 * Cshift(Foo, 0, -1); - ColourMatrix foo; - ColourMatrix bar; - ColourMatrix shifted1; - ColourMatrix shifted2; - ColourMatrix shifted3; - ColourMatrix foobar1; - ColourMatrix foobar2; - ColourMatrix mdiff,amdiff; - - peekSite(shifted1,Shifted,coor); - peekSite(shifted2,Foo,shiftcoor); - peekSite(shifted3,ShiftedCheck,coor); - peekSite(foo,Foo,coor); - - mdiff = shifted1-shifted2; - amdiff=adj(mdiff); - ColourMatrix prod = amdiff*mdiff; - Complex trprod = trace(prod); - Real Ttr=real(trprod); - double nn=Ttr; - if ( nn > 0 ) - cout<<"Shift real trace fail "< 0 ) - cout<<"Shift fail (shifted1/shifted2-ref) "< 0 ) - cout<<"Shift rb fail (shifted3/shifted2-ref) "<black + rShifted = Cshift(bFoo, dir, shift); // Shift black->red + ShiftedCheck = zero; + setCheckerboard(ShiftedCheck, bShifted); // Put them all together + setCheckerboard(ShiftedCheck, + rShifted); // and check the results (later) - /* - // Testing Smearing routine compilation, separate in a different file - GridCartesian Fine(latt_size,simd_layout,mpi_layout); - Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation - Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing); - SmearedConfiguration< PeriodicGimplR > SmartConf(&Fine, 3, StoutSmearing); - - std::cout< coor(4); + for (coor[3] = 0; coor[3] < latt_size[3] / mpi_layout[3]; coor[3]++) { + for (coor[2] = 0; coor[2] < latt_size[2] / mpi_layout[2]; + coor[2]++) { + for (coor[1] = 0; coor[1] < latt_size[1] / mpi_layout[1]; + coor[1]++) { + for (coor[0] = 0; coor[0] < latt_size[0] / mpi_layout[0]; + coor[0]++) { + std::complex diff; + + std::vector shiftcoor = coor; + shiftcoor[dir] = (shiftcoor[dir] + shift + latt_size[dir]) % + (latt_size[dir] / mpi_layout[dir]); + + std::vector rl(4); + for (int dd = 0; dd < 4; dd++) { + rl[dd] = latt_size[dd] / simd_layout[dd] / mpi_layout[dd]; + } + int lex = coor[0] % rl[0] + (coor[1] % rl[1]) * rl[0] + + (coor[2] % rl[2]) * rl[0] * rl[1] + + (coor[3] % rl[3]) * rl[0] * rl[1] * rl[2]; + lex += +1000 * (coor[0] / rl[0]) + + 1000 * (coor[1] / rl[1]) * simd_layout[0] + + 1000 * (coor[2] / rl[2]) * simd_layout[0] * + simd_layout[1] + + 1000 * (coor[3] / rl[3]) * simd_layout[0] * + simd_layout[1] * simd_layout[2]; + + int lex_coor = shiftcoor[0] % rl[0] + + (shiftcoor[1] % rl[1]) * rl[0] + + (shiftcoor[2] % rl[2]) * rl[0] * rl[1] + + (shiftcoor[3] % rl[3]) * rl[0] * rl[1] * rl[2]; + lex_coor += +1000 * (shiftcoor[0] / rl[0]) + + 1000 * (shiftcoor[1] / rl[1]) * simd_layout[0] + + 1000 * (shiftcoor[2] / rl[2]) * simd_layout[0] * + simd_layout[1] + + 1000 * (shiftcoor[3] / rl[3]) * simd_layout[0] * + simd_layout[1] * simd_layout[2]; + + ColourMatrix foo; + ColourMatrix bar; + ColourMatrix shifted1; + ColourMatrix shifted2; + ColourMatrix shifted3; + ColourMatrix foobar1; + ColourMatrix foobar2; + ColourMatrix mdiff, amdiff; + + peekSite(shifted1, Shifted, coor); + peekSite(shifted2, Foo, shiftcoor); + peekSite(shifted3, ShiftedCheck, coor); + peekSite(foo, Foo, coor); + + mdiff = shifted1 - shifted2; + amdiff = adj(mdiff); + ColourMatrix prod = amdiff * mdiff; + Complex trprod = trace(prod); + Real Ttr = real(trprod); + double nn = Ttr; + if (nn > 0) + cout << "Shift real trace fail " << coor[0] << coor[1] + << coor[2] << coor[3] << endl; + + for (int r = 0; r < 3; r++) { + for (int c = 0; c < 3; c++) { + diff = shifted1()()(r, c) - shifted2()()(r, c); + nn = real(conjugate(diff) * diff); + if (nn > 0) + cout << "Shift fail (shifted1/shifted2-ref) " << coor[0] + << coor[1] << coor[2] << coor[3] << " " + << shifted1()()(r, c) << " " << shifted2()()(r, c) + << " " << foo()()(r, c) << " lex expect " + << lex_coor << " lex " << lex << endl; + else if (0) + cout << "Shift pass 1vs2 " << coor[0] << coor[1] + << coor[2] << coor[3] << " " << shifted1()()(r, c) + << " " << shifted2()()(r, c) << " " + << foo()()(r, c) << " lex expect " << lex_coor + << " lex " << lex << endl; + } + } + + for (int r = 0; r < 3; r++) { + for (int c = 0; c < 3; c++) { + diff = shifted3()()(r, c) - shifted2()()(r, c); + nn = real(conjugate(diff) * diff); + if (nn > 0) + cout << "Shift rb fail (shifted3/shifted2-ref) " + << coor[0] << coor[1] << coor[2] << coor[3] << " " + << shifted3()()(r, c) << " " << shifted2()()(r, c) + << " " << foo()()(r, c) << " lex expect " + << lex_coor << " lex " << lex << endl; + else if (0) + cout << "Shift rb pass 3vs2 " << coor[0] << coor[1] + << coor[2] << coor[3] << " " << shifted3()()(r, c) + << " " << shifted2()()(r, c) << " " + << foo()()(r, c) << " lex expect " << lex_coor + << " lex " << lex << endl; + } + } + peekSite(bar, Bar, coor); + + peekSite(foobar1, FooBar, coor); + foobar2 = foo * bar; + for (int r = 0; r < Nc; r++) { + for (int c = 0; c < Nc; c++) { + diff = foobar2()()(r, c) - foobar1()()(r, c); + nrm = nrm + real(conjugate(diff) * diff); + } + } + } + } + } + } + if (Fine.IsBoss()) { + std::cout << GridLogMessage + << "LatticeColorMatrix * LatticeColorMatrix nrm diff = " + << nrm << std::endl; + } + } + } + + } // loop for lat + } // loop for omp + + /* + // Testing Smearing routine compilation, separate in a different file + GridCartesian Fine(latt_size,simd_layout,mpi_layout); + Smear_APE< PeriodicGimplR > APEsmearing; // periodic gauge implemetation + Smear_Stout< PeriodicGimplR > StoutSmearing(&APEsmearing); + SmearedConfiguration< PeriodicGimplR > SmartConf(&Fine, 3, StoutSmearing); + + std::cout<