diff --git a/benchmarks/Benchmark_dwf_ntpf b/benchmarks/Benchmark_dwf_ntpf deleted file mode 100755 index ab9999c9..00000000 Binary files a/benchmarks/Benchmark_dwf_ntpf and /dev/null differ diff --git a/benchmarks/Benchmark_zmm b/benchmarks/Benchmark_zmm deleted file mode 100755 index c7eebe18..00000000 Binary files a/benchmarks/Benchmark_zmm and /dev/null differ diff --git a/lib/qcd/hmc/integrators/Integrator.h b/lib/qcd/hmc/integrators/Integrator.h index db61b114..d604d67a 100644 --- a/lib/qcd/hmc/integrators/Integrator.h +++ b/lib/qcd/hmc/integrators/Integrator.h @@ -1,33 +1,34 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./lib/qcd/hmc/integrators/Integrator.h +Source file: ./lib/qcd/hmc/integrators/Integrator.h - 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 */ //-------------------------------------------------------------------- /*! @file Integrator.h * @brief Classes for the Molecular Dynamics integrator @@ -40,208 +41,201 @@ Author: paboyle #ifndef INTEGRATOR_INCLUDED #define INTEGRATOR_INCLUDED -//class Observer; +// class Observer; #include - namespace Grid{ - namespace QCD{ +namespace Grid { +namespace QCD { - struct IntegratorParameters{ +struct IntegratorParameters { + int Nexp; + int MDsteps; // number of outer steps + RealD trajL; // trajectory length + RealD stepsize; - int Nexp; - int MDsteps; // number of outer steps - RealD trajL; // trajectory length - RealD stepsize; - - IntegratorParameters(int MDsteps_, - RealD trajL_=1.0, - int Nexp_=12): - Nexp(Nexp_), - MDsteps(MDsteps_), - trajL(trajL_), - stepsize(trajL/MDsteps) - { - // empty body constructor - }; - - }; - - /*! @brief Class for Molecular Dynamics management */ - template - class Integrator { - - protected: - - typedef IntegratorParameters ParameterType; - - IntegratorParameters Params; - - const ActionSet as; - - int levels; // - double t_U; // Track time passing on each level and for U and for P - std::vector t_P; // - - GaugeField P; - - SmearingPolicy &Smearer; - - // Should match any legal (SU(n)) gauge field - // Need to use this template to match Ncol to pass to SU class - template void generate_momenta(Lattice< iVector< iScalar< iMatrix >, Nd> > & P,GridParallelRNG& pRNG){ - typedef Lattice< iScalar< iScalar< iMatrix > > > GaugeLinkField; - GaugeLinkField Pmu(P._grid); - Pmu = zero; - for(int mu=0;mu::GaussianLieAlgebraMatrix(pRNG, Pmu); - PokeIndex(P, Pmu, mu); - } - } - - - //ObserverList observers; // not yet - // typedef std::vector ObserverList; - // void register_observers(); - // void notify_observers(); - - void update_P(GaugeField&U, int level, double ep){ - t_P[level]+=ep; - update_P(P,U,level,ep); - - std::cout<is_smeared); - as[level].actions.at(a)->deriv(Us,force); // deriv should NOT include Ta - - std::cout<< GridLogIntegrator << "Smearing (on/off): "<is_smeared <is_smeared) Smearer.smeared_force(force); - force = Ta(force); - std::cout<< GridLogIntegrator << "Force average: "<< norm2(force)/(U._grid->gSites()) <(U, mu); - auto Pmu=PeekIndex(Mom, mu); - Umu = expMat(Pmu, ep, Params.Nexp)*Umu; - ProjectOnGroup(Umu); - PokeIndex(U, Umu, mu); - } - // Update the smeared fields, can be implemented as observer - Smearer.set_GaugeField(U); - } - - virtual void step (GaugeField& U,int level, int first,int last)=0; - -public: - - Integrator(GridBase* grid, - IntegratorParameters Par, - ActionSet & Aset, - SmearingPolicy &Sm): - Params(Par), - as(Aset), - P(grid), - levels(Aset.size()), - Smearer(Sm) - { - t_P.resize(levels,0.0); - t_U=0.0; - // initialization of smearer delegated outside of Integrator - }; - - virtual ~Integrator(){} - - //Initialization of momenta and actions - void refresh(GaugeField& U,GridParallelRNG &pRNG){ - std::cout<is_smeared); - as[level].actions.at(actionID)->refresh(Us, pRNG); - } - } - } - - // Calculate action - RealD S(GaugeField& U){// here also U not used - - LatticeComplex Hloc(U._grid); Hloc = zero; - // Momenta - for (int mu=0; mu (P, mu); - Hloc -= trace(Pmu*Pmu); - } - Complex Hsum = sum(Hloc); - - RealD H = Hsum.real(); - RealD Hterm; - std::cout<is_smeared); - Hterm = as[level].actions.at(actionID)->S(Us); - std::cout< +class Integrator { + protected: + typedef IntegratorParameters ParameterType; + + IntegratorParameters Params; + + const ActionSet as; + + int levels; // + double t_U; // Track time passing on each level and for U and for P + std::vector t_P; // + + GaugeField P; + + SmearingPolicy& Smearer; + + // Should match any legal (SU(n)) gauge field + // Need to use this template to match Ncol to pass to SU class + template + void generate_momenta(Lattice >, Nd> >& P, + GridParallelRNG& pRNG) { + typedef Lattice > > > GaugeLinkField; + GaugeLinkField Pmu(P._grid); + Pmu = zero; + for (int mu = 0; mu < Nd; mu++) { + SU::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); + PokeIndex(P, Pmu, mu); + } + } + + // ObserverList observers; // not yet + // typedef std::vector ObserverList; + // void register_observers(); + // void notify_observers(); + + void update_P(GaugeField& U, int level, double ep) { + t_P[level] += ep; + update_P(P, U, level, ep); + + std::cout << GridLogIntegrator << "[" << level << "] P " + << " dt " << ep << " : t_P " << t_P[level] << std::endl; + } + + void update_P(GaugeField& Mom, GaugeField& U, int level, double ep) { + // input U actually not used... + for (int a = 0; a < as[level].actions.size(); ++a) { + GaugeField force(U._grid); + GaugeField& Us = Smearer.get_U(as[level].actions.at(a)->is_smeared); + as[level].actions.at(a)->deriv(Us, force); // deriv should NOT include Ta + + std::cout << GridLogIntegrator + << "Smearing (on/off): " << as[level].actions.at(a)->is_smeared + << std::endl; + if (as[level].actions.at(a)->is_smeared) Smearer.smeared_force(force); + force = Ta(force); + std::cout << GridLogIntegrator + << "Force average: " << norm2(force) / (U._grid->gSites()) + << std::endl; + Mom -= force * ep; + } + } + + void update_U(GaugeField& U, double ep) { + update_U(P, U, ep); + + t_U += ep; + int fl = levels - 1; + std::cout << GridLogIntegrator << " " + << "[" << fl << "] U " + << " dt " << ep << " : t_U " << t_U << std::endl; + } + void update_U(GaugeField& Mom, GaugeField& U, double ep) { + // rewrite exponential to deal automatically with the lorentz index? + // GaugeLinkField Umu(U._grid); + // GaugeLinkField Pmu(U._grid); + for (int mu = 0; mu < Nd; mu++) { + auto Umu = PeekIndex(U, mu); + auto Pmu = PeekIndex(Mom, mu); + Umu = expMat(Pmu, ep, Params.Nexp) * Umu; + ProjectOnGroup(Umu); + PokeIndex(U, Umu, mu); + } + // Update the smeared fields, can be implemented as observer + Smearer.set_GaugeField(U); + } + + virtual void step(GaugeField& U, int level, int first, int last) = 0; + + public: + Integrator(GridBase* grid, IntegratorParameters Par, + ActionSet& Aset, SmearingPolicy& Sm) + : Params(Par), as(Aset), P(grid), levels(Aset.size()), Smearer(Sm) { + t_P.resize(levels, 0.0); + t_U = 0.0; + // initialization of smearer delegated outside of Integrator + }; + + virtual ~Integrator() {} + + // Initialization of momenta and actions + void refresh(GaugeField& U, GridParallelRNG& pRNG) { + std::cout << GridLogIntegrator << "Integrator refresh\n"; + generate_momenta(P, pRNG); + for (int level = 0; level < as.size(); ++level) { + for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) { + // get gauge field from the SmearingPolicy and + // based on the boolean is_smeared in actionID + GaugeField& Us = + Smearer.get_U(as[level].actions.at(actionID)->is_smeared); + as[level].actions.at(actionID)->refresh(Us, pRNG); + } + } + } + + // Calculate action + RealD S(GaugeField& U) { // here also U not used + + LatticeComplex Hloc(U._grid); + Hloc = zero; + // Momenta + for (int mu = 0; mu < Nd; mu++) { + auto Pmu = PeekIndex(P, mu); + Hloc -= trace(Pmu * Pmu); + } + Complex Hsum = sum(Hloc); + + RealD H = Hsum.real(); + RealD Hterm; + std::cout << GridLogMessage << "Momentum action H_p = " << H << "\n"; + + // Actions + for (int level = 0; level < as.size(); ++level) { + for (int actionID = 0; actionID < as[level].actions.size(); ++actionID) { + // get gauge field from the SmearingPolicy and + // based on the boolean is_smeared in actionID + GaugeField& Us = + Smearer.get_U(as[level].actions.at(actionID)->is_smeared); + Hterm = as[level].actions.at(actionID)->S(Us); + std::cout << GridLogMessage << "S Level " << level << " term " + << actionID << " H = " << Hterm << std::endl; + H += Hterm; + } + } + + return H; + } + + void integrate(GaugeField& U) { + // reset the clocks + t_U = 0; + for (int level = 0; level < as.size(); ++level) { + t_P[level] = 0; + } + + for (int step = 0; step < Params.MDsteps; ++step) { // MD step + int first_step = (step == 0); + int last_step = (step == Params.MDsteps - 1); + this->step(U, 0, first_step, last_step); + } + + // Check the clocks all match on all levels + for (int level = 0; level < as.size(); ++level) { + assert(fabs(t_U - t_P[level]) < 1.0e-6); // must be the same + std::cout << GridLogIntegrator << " times[" << level + << "]= " << t_P[level] << " " << t_U << std::endl; + } + + // and that we indeed got to the end of the trajectory + assert(fabs(t_U - Params.trajL) < 1.0e-6); + } +}; } } -#endif//INTEGRATOR_INCLUDED +#endif // INTEGRATOR_INCLUDED diff --git a/lib/qcd/utils/SUn.h b/lib/qcd/utils/SUn.h index e9403836..048e17c4 100644 --- a/lib/qcd/utils/SUn.h +++ b/lib/qcd/utils/SUn.h @@ -1,106 +1,124 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./lib/qcd/utils/SUn.h +Source file: ./lib/qcd/utils/SUn.h - 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 */ #ifndef QCD_UTIL_SUN_H #define QCD_UTIL_SUN_H namespace Grid { - namespace QCD { +namespace QCD { -template +template class SU { -public: + public: + static int generators(void) { return ncolour * ncolour - 1; } + static int su2subgroups(void) { return (ncolour * (ncolour - 1)) / 2; } - static int generators(void) { return ncolour*ncolour-1; } - static int su2subgroups(void) { return (ncolour*(ncolour-1))/2; } + template + using iSUnMatrix = iScalar > >; + template + using iSU2Matrix = iScalar > >; + template + using iSUnAdjointMatrix = iScalar > >; - template using iSUnMatrix = iScalar > > ; - template using iSU2Matrix = iScalar > > ; - ////////////////////////////////////////////////////////////////////////////////////////////////// - // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, SU<2>::LatticeMatrix etc... + // Types can be accessed as SU<2>::Matrix , SU<2>::vSUnMatrix, + // SU<2>::LatticeMatrix etc... ////////////////////////////////////////////////////////////////////////////////////////////////// - typedef iSUnMatrix Matrix; - typedef iSUnMatrix MatrixF; - typedef iSUnMatrix MatrixD; + typedef iSUnMatrix Matrix; + typedef iSUnMatrix MatrixF; + typedef iSUnMatrix MatrixD; - typedef iSUnMatrix vMatrix; - typedef iSUnMatrix vMatrixF; - typedef iSUnMatrix vMatrixD; + typedef iSUnMatrix vMatrix; + typedef iSUnMatrix vMatrixF; + typedef iSUnMatrix vMatrixD; - typedef Lattice LatticeMatrix; - typedef Lattice LatticeMatrixF; - typedef Lattice LatticeMatrixD; + // Actually the adjoint matrices are real... + // Consider this overhead... FIXME + typedef iSUnAdjointMatrix AMatrix; + typedef iSUnAdjointMatrix AMatrixF; + typedef iSUnAdjointMatrix AMatrixD; - typedef iSU2Matrix SU2Matrix; - typedef iSU2Matrix SU2MatrixF; - typedef iSU2Matrix SU2MatrixD; + typedef iSUnAdjointMatrix vAMatrix; + typedef iSUnAdjointMatrix vAMatrixF; + typedef iSUnAdjointMatrix vAMatrixD; - typedef iSU2Matrix vSU2Matrix; - typedef iSU2Matrix vSU2MatrixF; - typedef iSU2Matrix vSU2MatrixD; + typedef Lattice LatticeMatrix; + typedef Lattice LatticeMatrixF; + typedef Lattice LatticeMatrixD; - typedef Lattice LatticeSU2Matrix; - typedef Lattice LatticeSU2MatrixF; - typedef Lattice LatticeSU2MatrixD; + typedef iSU2Matrix SU2Matrix; + typedef iSU2Matrix SU2MatrixF; + typedef iSU2Matrix SU2MatrixD; + typedef iSU2Matrix vSU2Matrix; + typedef iSU2Matrix vSU2MatrixF; + typedef iSU2Matrix vSU2MatrixD; + + typedef Lattice LatticeSU2Matrix; + typedef Lattice LatticeSU2MatrixF; + typedef Lattice LatticeSU2MatrixD; //////////////////////////////////////////////////////////////////////// // There are N^2-1 generators for SU(N). // // We take a traceless hermitian generator basis as follows // - // * Normalisation: trace ta tb = 1/2 delta_ab - // + // * Normalisation: trace ta tb = 1/2 delta_ab = T_F delta_ab + // T_F = 1/2 for SU(N) groups + // // * Off diagonal // - pairs of rows i1,i2 behaving like pauli matrices signma_x, sigma_y - // - // - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] - // direct count off each row + // + // - there are (Nc-1-i1) slots for i2 on each row [ x 0 x ] + // direct count off each row // // - Sum of all pairs is Nc(Nc-1)/2: proof arithmetic series // - // (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 - // 1+ 2+ + + Nc-1 - // + // (Nc-1) + (Nc-2)+... 1 ==> Nc*(Nc-1)/2 + // 1+ 2+ + + Nc-1 + // // - There are 2 x Nc (Nc-1)/ 2 of these = Nc^2 - Nc // // - We enumerate the row-col pairs. - // - for each row col pair there is a (sigma_x) and a (sigma_y) like generator + // - for each row col pair there is a (sigma_x) and a (sigma_y) like + // generator // // - // t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => delta_{i,i1} delta_{j,i2} + delta_{i,i1} delta_{j,i2} - // t^a_ij = { in Nc(Nc-1)/2 ... Nc^(Nc-1) -1} => i delta_{i,i1} delta_{j,i2} - i delta_{i,i1} delta_{j,i2} - // + // t^a_ij = { in 0.. Nc(Nc-1)/2 -1} => 1/2(delta_{i,i1} delta_{j,i2} + + // delta_{i,i1} delta_{j,i2}) + // t^a_ij = { in Nc(Nc-1)/2 ... Nc(Nc-1) - 1} => i/2( delta_{i,i1} + // delta_{j,i2} - i delta_{i,i1} delta_{j,i2}) + // // * Diagonal; must be traceless and normalised - // - Sequence is + // - Sequence is // N (1,-1,0,0...) // N (1, 1,-2,0...) // N (1, 1, 1,-3,0...) @@ -109,558 +127,660 @@ public: // where 1/2 = N^2 (1+.. m^2)etc.... for the m-th diagonal generator // NB this gives the famous SU3 result for su2 index 8 // - // N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) + // N= sqrt(1/2 . 1/6 ) = 1/2 . 1/sqrt(3) // // ( 1 ) // ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 // ( -2) + // + // + // * Adjoint representation generators + // + // base for NxN hermitian traceless matrices + // normalized to 1: + // + // (e_Adj)^a = t^a / sqrt(T_F) + // + // then the real, antisymmetric generators for the adjoint representations + // are computed ( shortcut: e^a == (e_Adj)^a ) + // + // (iT_adj)^d_ab = i tr[e^a t^d e^b - t^d e^a e^b] + // //////////////////////////////////////////////////////////////////////// - template - static void generator(int lieIndex,iSUnMatrix &ta){ + template + static void generator(int lieIndex, iSUnMatrix &ta) { // map lie index to which type of generator int diagIndex; int su2Index; int sigxy; - int NNm1 = ncolour*(ncolour-1); - if ( lieIndex>= NNm1 ) { - diagIndex = lieIndex -NNm1; - generatorDiagonal(diagIndex,ta); + int NNm1 = ncolour * (ncolour - 1); + if (lieIndex >= NNm1) { + diagIndex = lieIndex - NNm1; + generatorDiagonal(diagIndex, ta); return; } - sigxy = lieIndex&0x1; - su2Index= lieIndex>>1; - if ( sigxy ) generatorSigmaY(su2Index,ta); - else generatorSigmaX(su2Index,ta); + sigxy = lieIndex & 0x1;//even or odd + su2Index = lieIndex >> 1; + if (sigxy) + generatorSigmaY(su2Index, ta); + else + generatorSigmaX(su2Index, ta); } - template - static void generatorSigmaX(int su2Index,iSUnMatrix &ta){ - ta=zero; - int i1,i2; - su2SubGroupIndex(i1,i2,su2Index); - ta()()(i1,i2)=1.0; - ta()()(i2,i1)=1.0; - ta= ta *0.5; + template + static void generatorSigmaX(int su2Index, iSUnMatrix &ta) { + ta = zero; + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = 1.0; + ta()()(i2, i1) = 1.0; + ta = ta * 0.5; } - template - static void generatorSigmaY(int su2Index,iSUnMatrix &ta){ - ta=zero; - cplx i(0.0,1.0); - int i1,i2; - su2SubGroupIndex(i1,i2,su2Index); - ta()()(i1,i2)=-i; - ta()()(i2,i1)= i; - ta= ta *0.5; + template + static void generatorSigmaY(int su2Index, iSUnMatrix &ta) { + ta = zero; + cplx i(0.0, 1.0); + int i1, i2; + su2SubGroupIndex(i1, i2, su2Index); + ta()()(i1, i2) = -i; + ta()()(i2, i1) = i; + ta = ta * 0.5; } - template - static void generatorDiagonal(int diagIndex,iSUnMatrix &ta){ - ta=zero; - int trsq=0; - int last=diagIndex+1; - for(int i=0;i<=diagIndex;i++){ - ta()()(i,i) = 1.0; - trsq++; + template + static void generatorDiagonal(int diagIndex, iSUnMatrix &ta) { + // diag ({1, 1, ..., 1}(k-times), -k, 0, 0, ...) + ta = zero; + int k = diagIndex + 1;// diagIndex starts from 0 + for (int i = 0; i <= diagIndex; i++) {// k iterations + ta()()(i, i) = 1.0; } - ta()()(last,last) = -last; - trsq+=last*last; - RealD nrm = 1.0/std::sqrt(2.0*trsq); - ta = ta *nrm; + ta()()(k,k) = -k;//indexing starts from 0 + RealD nrm = 1.0 / std::sqrt(2.0 * k*(k+1)); + ta = ta * nrm; } + template + static void generatorAdjoint(int Index, iSUnAdjointMatrix &iAdjTa){ + // returns i(T_Adj)^index necessary for the projectors + // see definitions above + iAdjTa = zero; + Vector< iSUnMatrix > ta(ncolour*ncolour -1); + iSUnMatrix tmp; + + // FIXME not very efficient to get all the generators everytime + for (int a = 0; a < (ncolour * ncolour - 1); a++) + generator(a, ta[a]); + + + for (int a = 0; a < (ncolour*ncolour - 1); a++){ + tmp = ta[a] * ta[Index] - ta[Index] * ta[a]; + for (int b = 0; b < (ncolour*ncolour - 1); b++){ + iSUnMatrix tmp1 = 2.0 * tmp * ta[b]; // 2.0 from the normalization + Complex iTr = TensorRemove(timesI(trace(tmp1))); + iAdjTa()()(b,a) = iTr; + } + } + + + } //////////////////////////////////////////////////////////////////////// // Map a su2 subgroup number to the pair of rows that are non zero //////////////////////////////////////////////////////////////////////// - static void su2SubGroupIndex(int &i1,int &i2,int su2_index){ + static void su2SubGroupIndex(int &i1, int &i2, int su2_index) { + assert((su2_index >= 0) && (su2_index < (ncolour * (ncolour - 1)) / 2)); - assert( (su2_index>=0) && (su2_index< (ncolour*(ncolour-1))/2) ); - - int spare=su2_index; - for(i1=0;spare >= (ncolour-1-i1);i1++ ){ - spare = spare - (ncolour-1-i1); // remove the Nc-1-i1 terms + int spare = su2_index; + for (i1 = 0; spare >= (ncolour - 1 - i1); i1++) { + spare = spare - (ncolour - 1 - i1); // remove the Nc-1-i1 terms } - i2=i1+1+spare; + i2 = i1 + 1 + spare; } ////////////////////////////////////////////////////////////////////////////////////////// // Pull out a subgroup and project on to real coeffs x pauli basis ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Extract( Lattice > &Determinant, - Lattice > &subgroup, - const Lattice > &source, - int su2_index) - { + template + static void su2Extract(Lattice > &Determinant, + Lattice > &subgroup, + const Lattice > &source, + int su2_index) { GridBase *grid(source._grid); - conformable(subgroup,source); - conformable(subgroup,Determinant); - int i0,i1; - su2SubGroupIndex(i0,i1,su2_index); + conformable(subgroup, source); + conformable(subgroup, Determinant); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - subgroup._odata[ss]()()(0,0) = source._odata[ss]()()(i0,i0); - subgroup._odata[ss]()()(0,1) = source._odata[ss]()()(i0,i1); - subgroup._odata[ss]()()(1,0) = source._odata[ss]()()(i1,i0); - subgroup._odata[ss]()()(1,1) = source._odata[ss]()()(i1,i1); + PARALLEL_FOR_LOOP + for (int ss = 0; ss < grid->oSites(); ss++) { + subgroup._odata[ss]()()(0, 0) = source._odata[ss]()()(i0, i0); + subgroup._odata[ss]()()(0, 1) = source._odata[ss]()()(i0, i1); + subgroup._odata[ss]()()(1, 0) = source._odata[ss]()()(i1, i0); + subgroup._odata[ss]()()(1, 1) = source._odata[ss]()()(i1, i1); iSU2Matrix Sigma = subgroup._odata[ss]; - Sigma = Sigma-adj(Sigma)+trace(adj(Sigma)); + Sigma = Sigma - adj(Sigma) + trace(adj(Sigma)); subgroup._odata[ss] = Sigma; // this should be purely real - Determinant._odata[ss] = Sigma()()(0,0)*Sigma()()(1,1) - - Sigma()()(0,1)*Sigma()()(1,0); - + Determinant._odata[ss] = + Sigma()()(0, 0) * Sigma()()(1, 1) - Sigma()()(0, 1) * Sigma()()(1, 0); } - } ////////////////////////////////////////////////////////////////////////////////////////// // Set matrix to one and insert a pauli subgroup ////////////////////////////////////////////////////////////////////////////////////////// - template - static void su2Insert( const Lattice > &subgroup, - Lattice > &dest, - int su2_index) - { + template + static void su2Insert(const Lattice > &subgroup, + Lattice > &dest, int su2_index) { GridBase *grid(dest._grid); - conformable(subgroup,dest); - int i0,i1; - su2SubGroupIndex(i0,i1,su2_index); + conformable(subgroup, dest); + int i0, i1; + su2SubGroupIndex(i0, i1, su2_index); - dest = 1.0; // start out with identity -PARALLEL_FOR_LOOP - for(int ss=0;ssoSites();ss++){ - dest._odata[ss]()()(i0,i0) = subgroup._odata[ss]()()(0,0); - dest._odata[ss]()()(i0,i1) = subgroup._odata[ss]()()(0,1); - dest._odata[ss]()()(i1,i0) = subgroup._odata[ss]()()(1,0); - dest._odata[ss]()()(i1,i1) = subgroup._odata[ss]()()(1,1); + dest = 1.0; // start out with identity + PARALLEL_FOR_LOOP + for (int ss = 0; ss < grid->oSites(); ss++) { + dest._odata[ss]()()(i0, i0) = subgroup._odata[ss]()()(0, 0); + dest._odata[ss]()()(i0, i1) = subgroup._odata[ss]()()(0, 1); + dest._odata[ss]()()(i1, i0) = subgroup._odata[ss]()()(1, 0); + dest._odata[ss]()()(i1, i1) = subgroup._odata[ss]()()(1, 1); } } - /////////////////////////////////////////////// - // Generate e^{ Re Tr Staple Link} dlink + // Generate e^{ Re Tr Staple Link} dlink // - // *** Note Staple should be appropriate linear compbination between all staples. + // *** Note Staple should be appropriate linear compbination between all + // staples. // *** If already by beta pass coefficient 1.0. - // *** This routine applies the additional 1/Nc factor that comes after trace in action. + // *** This routine applies the additional 1/Nc factor that comes after trace + // in action. // /////////////////////////////////////////////// - static void SubGroupHeatBath( GridSerialRNG &sRNG, - GridParallelRNG &pRNG, - RealD beta,//coeff multiplying staple in action (with no 1/Nc) - LatticeMatrix &link, - const LatticeMatrix &barestaple, // multiplied by action coeffs so th - int su2_subgroup, - int nheatbath, - LatticeInteger &wheremask) - { + static void SubGroupHeatBath( + GridSerialRNG &sRNG, GridParallelRNG &pRNG, + RealD beta, // coeff multiplying staple in action (with no 1/Nc) + LatticeMatrix &link, + const LatticeMatrix &barestaple, // multiplied by action coeffs so th + int su2_subgroup, int nheatbath, LatticeInteger &wheremask) { GridBase *grid = link._grid; - int ntrials=0; - int nfails=0; - const RealD twopi=2.0*M_PI; + int ntrials = 0; + int nfails = 0; + const RealD twopi = 2.0 * M_PI; - LatticeMatrix staple(grid); + LatticeMatrix staple(grid); - staple = barestaple * (beta/ncolour); + staple = barestaple * (beta / ncolour); - LatticeMatrix V(grid); V = link*staple; + LatticeMatrix V(grid); + V = link * staple; // Subgroup manipulation in the lie algebra space - LatticeSU2Matrix u(grid); // Kennedy pendleton "u" real projected normalised Sigma + LatticeSU2Matrix u( + grid); // Kennedy pendleton "u" real projected normalised Sigma LatticeSU2Matrix uinv(grid); - LatticeSU2Matrix ua(grid); // a in pauli form - LatticeSU2Matrix b(grid); // rotated matrix after hb + LatticeSU2Matrix ua(grid); // a in pauli form + LatticeSU2Matrix b(grid); // rotated matrix after hb // Some handy constant fields - LatticeComplex ones (grid); ones = 1.0; - LatticeComplex zeros(grid); zeros=zero; - LatticeReal rones (grid); rones = 1.0; - LatticeReal rzeros(grid); rzeros=zero; - LatticeComplex udet(grid); // determinant of real(staple) - LatticeInteger mask_true (grid); mask_true = 1; - LatticeInteger mask_false(grid); mask_false= 0; + LatticeComplex ones(grid); + ones = 1.0; + LatticeComplex zeros(grid); + zeros = zero; + LatticeReal rones(grid); + rones = 1.0; + LatticeReal rzeros(grid); + rzeros = zero; + LatticeComplex udet(grid); // determinant of real(staple) + LatticeInteger mask_true(grid); + mask_true = 1; + LatticeInteger mask_false(grid); + mask_false = 0; - /* -PLB 156 P393 (1985) (Kennedy and Pendleton) + /* + PLB 156 P393 (1985) (Kennedy and Pendleton) -Note: absorb "beta" into the def of sigma compared to KP paper; staple -passed to this routine has "beta" already multiplied in + Note: absorb "beta" into the def of sigma compared to KP paper; staple + passed to this routine has "beta" already multiplied in -Action linear in links h and of form: + Action linear in links h and of form: - beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) + beta S = beta Sum_p (1 - 1/Nc Re Tr Plaq ) -Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " + Writing Sigma = 1/Nc (beta Sigma') where sum over staples is "Sigma' " - beta S = const - beta/Nc Re Tr h Sigma' - = const - Re Tr h Sigma - -Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex arbitrary. + beta S = const - beta/Nc Re Tr h Sigma' + = const - Re Tr h Sigma - Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij - Re Tr h Sigma = 2 h_j Re Sigma_j + Decompose h and Sigma into (1, sigma_j) ; h_i real, h^2=1, Sigma_i complex + arbitrary. -Normalised re Sigma_j = xi u_j + Tr h Sigma = h_i Sigma_j Tr (sigma_i sigma_j) = h_i Sigma_j 2 delta_ij + Re Tr h Sigma = 2 h_j Re Sigma_j -With u_j a unit vector and U can be in SU(2); + Normalised re Sigma_j = xi u_j -Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) + With u_j a unit vector and U can be in SU(2); -4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + Re Tr h Sigma = 2 h_j Re Sigma_j = 2 xi (h.u) - xi = sqrt(Det)/2; + 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] -Write a= u h in SU(2); a has pauli decomp a_j; + xi = sqrt(Det)/2; -Note: Product b' xi is unvariant because scaling Sigma leaves - normalised vector "u" fixed; Can rescale Sigma so b' = 1. - */ + Write a= u h in SU(2); a has pauli decomp a_j; + + Note: Product b' xi is unvariant because scaling Sigma leaves + normalised vector "u" fixed; Can rescale Sigma so b' = 1. + */ //////////////////////////////////////////////////////// // Real part of Pauli decomposition // Note a subgroup can project to zero in cold start //////////////////////////////////////////////////////// - su2Extract(udet,u,V,su2_subgroup); + su2Extract(udet, u, V, su2_subgroup); ////////////////////////////////////////////////////// // Normalising this vector if possible; else identity ////////////////////////////////////////////////////// - LatticeComplex xi(grid); + LatticeComplex xi(grid); LatticeSU2Matrix lident(grid); - - SU2Matrix ident = Complex(1.0); - SU2Matrix pauli1; SU<2>::generator(0,pauli1); - SU2Matrix pauli2; SU<2>::generator(1,pauli2); - SU2Matrix pauli3; SU<2>::generator(2,pauli3); - pauli1 = timesI(pauli1)*2.0; - pauli2 = timesI(pauli2)*2.0; - pauli3 = timesI(pauli3)*2.0; + + SU2Matrix ident = Complex(1.0); + SU2Matrix pauli1; + SU<2>::generator(0, pauli1); + SU2Matrix pauli2; + SU<2>::generator(1, pauli2); + SU2Matrix pauli3; + SU<2>::generator(2, pauli3); + pauli1 = timesI(pauli1) * 2.0; + pauli2 = timesI(pauli2) * 2.0; + pauli3 = timesI(pauli3) * 2.0; LatticeComplex cone(grid); LatticeReal adet(grid); adet = abs(toReal(udet)); - lident=Complex(1.0); - cone =Complex(1.0); - Real machine_epsilon=1.0e-7; - u = where(adet>machine_epsilon,u,lident); - udet= where(adet>machine_epsilon,udet,cone); + lident = Complex(1.0); + cone = Complex(1.0); + Real machine_epsilon = 1.0e-7; + u = where(adet > machine_epsilon, u, lident); + udet = where(adet > machine_epsilon, udet, cone); - xi = 0.5*sqrt(udet); //4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] - u = 0.5*u*pow(xi,-1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] + xi = 0.5 * sqrt(udet); // 4xi^2 = Det [ Sig - Sig^dag + 1 Tr Sigdag] + u = 0.5 * u * + pow(xi, -1.0); // u = 1/2xi [ Sig - Sig^dag + 1 Tr Sigdag] // Debug test for sanity - uinv=adj(u); - b=u*uinv-1.0; - assert(norm2(b)<1.0e-4); + uinv = adj(u); + b = u * uinv - 1.0; + assert(norm2(b) < 1.0e-4); - /* -Measure: Haar measure dh has d^4a delta(1-|a^2|) -In polars: - da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) - = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + r) ) - = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) + /* + Measure: Haar measure dh has d^4a delta(1-|a^2|) + In polars: + da = da0 r^2 sin theta dr dtheta dphi delta( 1 - r^2 -a0^2) + = da0 r^2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r)(sqrt(1-a0^) + + r) ) + = da0 r/2 sin theta dr dtheta dphi delta( (sqrt(1-a0^) - r) ) -Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters through xi - = e^{2 xi (h.u)} dh - = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi h2u2}.e^{2 xi h3u3} dh + Action factor Q(h) dh = e^-S[h] dh = e^{ xi Tr uh} dh // beta enters + through xi + = e^{2 xi (h.u)} dh + = e^{2 xi h0u0}.e^{2 xi h1u1}.e^{2 xi + h2u2}.e^{2 xi h3u3} dh -Therefore for each site, take xi for that site -i) generate |a0|<1 with dist - (1-a0^2)^0.5 e^{2 xi a0 } da0 + Therefore for each site, take xi for that site + i) generate |a0|<1 with dist + (1-a0^2)^0.5 e^{2 xi a0 } da0 -Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc factor in Chroma ] -A. Generate two uniformly distributed pseudo-random numbers R and R', R'', R''' in the unit interval; -B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha; -C. Set C = cos^2(2pi R"), with R" another uniform random number in [0,1] ; -D. Set A = XC; -E. Let d = X'+A; -F. If R'''^2 :> 1 - 0.5 d, go back to A; -G. Set a0 = 1 - d; + Take alpha = 2 xi = 2 xi [ recall 2 beta/Nc unmod staple norm]; hence 2.0/Nc + factor in Chroma ] + A. Generate two uniformly distributed pseudo-random numbers R and R', R'', + R''' in the unit interval; + B. Set X = -(ln R)/alpha, X' =-(ln R')/alpha; + C. Set C = cos^2(2pi R"), with R" another uniform random number in [0,1] ; + D. Set A = XC; + E. Let d = X'+A; + F. If R'''^2 :> 1 - 0.5 d, go back to A; + G. Set a0 = 1 - d; -Note that in step D setting B ~ X - A and using B in place of A in step E will generate a second independent a 0 value. - */ + Note that in step D setting B ~ X - A and using B in place of A in step E will + generate a second independent a 0 value. + */ ///////////////////////////////////////////////////////// // count the number of sites by picking "1"'s out of hat ///////////////////////////////////////////////////////// - Integer hit=0; + Integer hit = 0; LatticeReal rtmp(grid); - rtmp=where(wheremask,rones,rzeros); + rtmp = where(wheremask, rones, rzeros); RealD numSites = sum(rtmp); RealD numAccepted; - LatticeInteger Accepted(grid); Accepted = zero; + LatticeInteger Accepted(grid); + Accepted = zero; LatticeInteger newlyAccepted(grid); - - std::vector xr(4,grid); - std::vector a(4,grid); - LatticeReal d(grid); d=zero; + + std::vector xr(4, grid); + std::vector a(4, grid); + LatticeReal d(grid); + d = zero; LatticeReal alpha(grid); // std::cout< 1 - 0.5 d, go back to A; - LatticeReal thresh(grid); thresh = 1.0-d*0.5; - xrsq = xr[0]*xr[0]; - LatticeInteger ione(grid); ione = 1; - LatticeInteger izero(grid); izero=zero; + d = where(Accepted, d, xr[2] + xr[1] * xr[3]); - newlyAccepted = where ( xrsq < thresh,ione,izero); - Accepted = where ( newlyAccepted, newlyAccepted,Accepted); - Accepted = where ( wheremask, Accepted,izero); + // F. If R'''^2 :> 1 - 0.5 d, go back to A; + LatticeReal thresh(grid); + thresh = 1.0 - d * 0.5; + xrsq = xr[0] * xr[0]; + LatticeInteger ione(grid); + ione = 1; + LatticeInteger izero(grid); + izero = zero; + + newlyAccepted = where(xrsq < thresh, ione, izero); + Accepted = where(newlyAccepted, newlyAccepted, Accepted); + Accepted = where(wheremask, Accepted, izero); // FIXME need an iSum for integer to avoid overload on return type?? - rtmp=where(Accepted,rones,rzeros); + rtmp = where(Accepted, rones, rzeros); numAccepted = sum(rtmp); hit++; - } while ( (numAccepted < numSites) && ( hit < nheatbath) ); + } while ((numAccepted < numSites) && (hit < nheatbath)); // G. Set a0 = 1 - d; - a[0]=zero; - a[0]=where(wheremask,1.0-d,a[0]); + a[0] = zero; + a[0] = where(wheremask, 1.0 - d, a[0]); ////////////////////////////////////////// // ii) generate a_i uniform on two sphere radius (1-a0^2)^0.5 ////////////////////////////////////////// LatticeReal a123mag(grid); - a123mag = sqrt(abs(1.0-a[0]*a[0])); + a123mag = sqrt(abs(1.0 - a[0] * a[0])); - LatticeReal cos_theta (grid); - LatticeReal sin_theta (grid); - LatticeReal phi (grid); + LatticeReal cos_theta(grid); + LatticeReal sin_theta(grid); + LatticeReal phi(grid); - random(pRNG,phi); phi = phi * twopi; // uniform in [0,2pi] - random(pRNG,cos_theta); cos_theta=(cos_theta*2.0)-1.0; // uniform in [-1,1] - sin_theta = sqrt(abs(1.0-cos_theta*cos_theta)); + random(pRNG, phi); + phi = phi * twopi; // uniform in [0,2pi] + random(pRNG, cos_theta); + cos_theta = (cos_theta * 2.0) - 1.0; // uniform in [-1,1] + sin_theta = sqrt(abs(1.0 - cos_theta * cos_theta)); - a[1] = a123mag * sin_theta * cos(phi); - a[2] = a123mag * sin_theta * sin(phi); - a[3] = a123mag * cos_theta; - - ua = toComplex(a[0])*ident - + toComplex(a[1])*pauli1 - + toComplex(a[2])*pauli2 - + toComplex(a[3])*pauli3; + a[1] = a123mag * sin_theta * cos(phi); + a[2] = a123mag * sin_theta * sin(phi); + a[3] = a123mag * cos_theta; + + ua = toComplex(a[0]) * ident + toComplex(a[1]) * pauli1 + + toComplex(a[2]) * pauli2 + toComplex(a[3]) * pauli3; b = 1.0; - b = where(wheremask,uinv*ua,b); - su2Insert(b,V,su2_subgroup); + b = where(wheremask, uinv * ua, b); + su2Insert(b, V, su2_subgroup); - //mask the assignment back based on Accptance - link = where(Accepted,V * link,link); + // mask the assignment back based on Accptance + link = where(Accepted, V * link, link); ////////////////////////////// // Debug Checks // SU2 check - LatticeSU2Matrix check(grid); // rotated matrix after hb + LatticeSU2Matrix check(grid); // rotated matrix after hb u = zero; check = ua * adj(ua) - 1.0; - check = where(Accepted,check,u); - assert(norm2(check)<1.0e-4); + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); + + check = b * adj(b) - 1.0; + check = where(Accepted, check, u); + assert(norm2(check) < 1.0e-4); - check = b*adj(b)-1.0; - check = where(Accepted,check,u); - assert(norm2(check)<1.0e-4); - LatticeMatrix Vcheck(grid); Vcheck = zero; - Vcheck = where(Accepted,V*adj(V) - 1.0,Vcheck); + Vcheck = where(Accepted, V * adj(V) - 1.0, Vcheck); // std::cout< - static void LieRandomize(GridParallelRNG &pRNG,LatticeMatrixType &out,double scale=1.0){ + template + static void LieRandomize(GridParallelRNG &pRNG, LatticeMatrixType &out, + double scale = 1.0) { GridBase *grid = out._grid; - + typedef typename LatticeMatrixType::vector_type vector_type; typedef typename LatticeMatrixType::scalar_type scalar_type; - + typedef iSinglet vTComplexType; - + typedef Lattice LatticeComplexType; - typedef typename GridTypeMapper::scalar_object MatrixType; - - LatticeComplexType ca (grid); - LatticeMatrixType lie(grid); - LatticeMatrixType la (grid); - ComplexD ci(0.0,scale); - ComplexD cone(1.0,0.0); + typedef typename GridTypeMapper< + typename LatticeMatrixType::vector_object>::scalar_object MatrixType; + + LatticeComplexType ca(grid); + LatticeMatrixType lie(grid); + LatticeMatrixType la(grid); + ComplexD ci(0.0, scale); + ComplexD cone(1.0, 0.0); MatrixType ta; - lie=zero; - for(int a=0;a - static void HotConfiguration(GridParallelRNG &pRNG,GaugeField &out){ + static void FundamentalLieAlgebraMatrix(Vector &h, LatticeMatrix &out, + Real scale = 1.0) { + GridBase *grid = out._grid; + LatticeMatrix la(grid); + Matrix ta; + + out = zero; + for (int a = 0; a < generators(); a++) { + generator(a, ta); + la = Complex(0.0, h[a]) * scale * ta; + out += la; + } + } + + template + static void HotConfiguration(GridParallelRNG &pRNG, GaugeField &out) { typedef typename GaugeField::vector_type vector_type; typedef iSUnMatrix vMatrixType; typedef Lattice LatticeMatrixType; - + LatticeMatrixType Umu(out._grid); - for(int mu=0;mu(out,Umu,mu); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 1.0); + PokeIndex(out, Umu, mu); } } - static void TepidConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){ + static void TepidConfiguration(GridParallelRNG &pRNG, + LatticeGaugeField &out) { LatticeMatrix Umu(out._grid); - for(int mu=0;mu(out,Umu,mu); + for (int mu = 0; mu < Nd; mu++) { + LieRandomize(pRNG, Umu, 0.01); + PokeIndex(out, Umu, mu); } } - static void ColdConfiguration(GridParallelRNG &pRNG,LatticeGaugeField &out){ + static void ColdConfiguration(GridParallelRNG &pRNG, LatticeGaugeField &out) { LatticeMatrix Umu(out._grid); - Umu=1.0; - for(int mu=0;mu(out,Umu,mu); + Umu = 1.0; + for (int mu = 0; mu < Nd; mu++) { + PokeIndex(out, Umu, mu); } } - static void taProj( const LatticeMatrix &in, LatticeMatrix &out){ + static void taProj(const LatticeMatrix &in, LatticeMatrix &out) { out = Ta(in); } - template - static void taExp( const LatticeMatrixType &x, LatticeMatrixType &ex){ - typedef typename LatticeMatrixType::scalar_type ComplexType; - + template + static void taExp(const LatticeMatrixType &x, LatticeMatrixType &ex) { + typedef typename LatticeMatrixType::scalar_type ComplexType; + LatticeMatrixType xn(x._grid); RealD nfac = 1.0; xn = x; - ex =xn+ComplexType(1.0); // 1+x + ex = xn + ComplexType(1.0); // 1+x // Do a 12th order exponentiation - for(int i=2; i <= 12; ++i) - { - nfac = nfac/RealD(i); //1/2, 1/2.3 ... - xn = xn * x; // x2, x3,x4.... - ex = ex+ xn*nfac;// x2/2!, x3/3!.... + for (int i = 2; i <= 12; ++i) { + nfac = nfac / RealD(i); // 1/2, 1/2.3 ... + xn = xn * x; // x2, x3,x4.... + ex = ex + xn * nfac; // x2/2!, x3/3!.... } } - }; - typedef SU<2> SU2; - typedef SU<3> SU3; - typedef SU<4> SU4; - typedef SU<5> SU5; - - } +typedef SU<2> SU2; +typedef SU<3> SU3; +typedef SU<4> SU4; +typedef SU<5> SU5; +} } #endif diff --git a/tests/Test_contfrac_force.cc b/tests/Test_contfrac_force.cc index b8a26d4d..0695537c 100644 --- a/tests/Test_contfrac_force.cc +++ b/tests/Test_contfrac_force.cc @@ -96,7 +96,7 @@ int main (int argc, char ** argv) for(int mu=0;mu(mom,mommu,mu); diff --git a/tests/Test_dwf_force.cc b/tests/Test_dwf_force.cc index 389a7318..734d03a4 100644 --- a/tests/Test_dwf_force.cc +++ b/tests/Test_dwf_force.cc @@ -96,7 +96,7 @@ int main (int argc, char ** argv) for(int mu=0;mu(mom,mommu,mu); diff --git a/tests/Test_dwf_gpforce.cc b/tests/Test_dwf_gpforce.cc index 4a8ceb5b..e336a924 100644 --- a/tests/Test_dwf_gpforce.cc +++ b/tests/Test_dwf_gpforce.cc @@ -113,7 +113,7 @@ int main (int argc, char ** argv) for(int mu=0;mu(mom,mommu,mu); diff --git a/tests/Test_gpdwf_force.cc b/tests/Test_gpdwf_force.cc index c600a2ff..3df2dc36 100644 --- a/tests/Test_gpdwf_force.cc +++ b/tests/Test_gpdwf_force.cc @@ -100,7 +100,7 @@ int main (int argc, char ** argv) for(int mu=0;mu(mom,mommu,mu); @@ -169,7 +169,7 @@ int main (int argc, char ** argv) // // Pmu = zero; // for(int mu=0;mu::GaussianLieAlgebraMatrix(pRNG, Pmu); + // SU::GaussianFundamentalLieAlgebraMatrix(pRNG, Pmu); // PokeIndex(P, Pmu, mu); // } // diff --git a/tests/Test_lie_generators.cc b/tests/Test_lie_generators.cc index 02620718..11dba39c 100644 --- a/tests/Test_lie_generators.cc +++ b/tests/Test_lie_generators.cc @@ -1,70 +1,74 @@ - /************************************************************************************* +/************************************************************************************* - Grid physics library, www.github.com/paboyle/Grid +Grid physics library, www.github.com/paboyle/Grid - Source file: ./tests/Test_lie_generators.cc +Source file: ./tests/Test_lie_generators.cc - Copyright (C) 2015 +Copyright (C) 2015 Author: Azusa Yamaguchi 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 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 #include -#include #include +#include using namespace std; using namespace Grid; using namespace Grid::QCD; +int main(int argc, char** argv) { + Grid_init(&argc, &argv); -int main (int argc, char ** argv) -{ - Grid_init(&argc,&argv); + std::vector latt({4, 4, 4, 8}); + GridCartesian* grid = SpaceTimeGrid::makeFourDimGrid( + latt, GridDefaultSimd(Nd, vComplex::Nsimd()), GridDefaultMpi()); - std::vector latt({4,4,4,8}); - GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, - GridDefaultSimd(Nd,vComplex::Nsimd()), - GridDefaultMpi()); - - GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); + GridRedBlackCartesian* rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); - std::cout<(mom,mommu,mu); diff --git a/tests/Test_rect_force.cc b/tests/Test_rect_force.cc index f437fa21..32a8a7cd 100644 --- a/tests/Test_rect_force.cc +++ b/tests/Test_rect_force.cc @@ -81,7 +81,7 @@ int main (int argc, char ** argv) for(int mu=0;mu(mom,mommu,mu); diff --git a/tests/Test_wilson_force.cc b/tests/Test_wilson_force.cc index 751c0896..7ac21852 100644 --- a/tests/Test_wilson_force.cc +++ b/tests/Test_wilson_force.cc @@ -93,7 +93,7 @@ int main (int argc, char ** argv) for(int mu=0;mu