From 21f214d6c90fe092d44a1d757e7949a1919cf0c1 Mon Sep 17 00:00:00 2001 From: Azusa Yamaguchi Date: Sun, 14 Jun 2015 00:50:59 +0100 Subject: [PATCH] Apply a heatbath sweep --- tests/Test_lie_generators.cc | 271 +++++++++-------------------------- 1 file changed, 68 insertions(+), 203 deletions(-) diff --git a/tests/Test_lie_generators.cc b/tests/Test_lie_generators.cc index b29682a0..91f2d33f 100644 --- a/tests/Test_lie_generators.cc +++ b/tests/Test_lie_generators.cc @@ -2,226 +2,91 @@ #include #include +#include using namespace std; using namespace Grid; using namespace Grid::QCD; -class suN { -public: - - static int generators(int ncolour) { return ncolour*ncolour-1; } - static int su2subgroups(int ncolour) { return (ncolour*(ncolour-1))/2; } - - template using suNmatrix = iScalar > > ; - - //////////////////////////////////////////////////////////////////////// - // 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 - // - // * 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 - // - // - 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 - // - // - 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 - // - // - // 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} - // - // * Diagonal; must be traceless and normalised - // - Sequence is - // N (1,-1,0,0...) - // N (1, 1,-2,0...) - // N (1, 1, 1,-3,0...) - // N (1, 1, 1, 1,-4,0...) - // - // 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) - // - // ( 1 ) - // ( 1 ) / sqrt(3) /2 = 1/2 lambda_8 - // ( -2) - //////////////////////////////////////////////////////////////////////// - template - static void suNgenerator(int lieIndex,suNmatrix &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; - suNgeneratorDiagonal(diagIndex,ta); - return; - } - sigxy = lieIndex&0x1; - su2Index= lieIndex>>1; - if ( sigxy ) suNgeneratorSigmaY(su2Index,ta); - else suNgeneratorSigmaX(su2Index,ta); - } - template - static void suNgeneratorSigmaX(int su2Index,suNmatrix &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 suNgeneratorSigmaY(int su2Index,suNmatrix &ta){ - ta=zero; - Complex 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 suNgeneratorDiagonal(int diagIndex,suNmatrix &ta){ - ta=zero; - int trsq=0; - int last=diagIndex+1; - for(int i=0;i<=diagIndex;i++){ - ta()()(i,i) = 1.0; - trsq++; - } - ta()()(last,last) = -last; - trsq+=last*last; - RealD nrm = 1.0/std::sqrt(2.0*trsq); - ta = ta *nrm; - } - //////////////////////////////////////////////////////////////////////// - // Map a - // - //////////////////////////////////////////////////////////////////////// - template - static void su2SubGroupIndex(int &i1,int &i2,int su2_index){ - - 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 - } - i2=i1+1+spare; - } - template - static void su2Extract(std::vector &r,const Lattice > &source, int su2_index) - { - assert(r.size() == 4); // store in 4 real parts - - for(int i=0;i<4;i++){ - conformable(r[i],source); - } - - int i1,i2; - su2SubGroupIndex(i1,i2,su2_index); - - /* Compute the b(k) of A_SU(2) = b0 + i sum_k bk sigma_k */ - r[0] = real(source()()(i1,i1) + source()()(i2,i2)); - r[1] = imag(source()()(i1,i2) + source()()(i2,i1)); - r[2] = real(source()()(i1,i2) - source()()(i2,i1)); - r[3] = imag(source()()(i1,i1) - source()()(i2,i2)); - } - - - template static void printGenerators(void) - { - for(int gen=0;gen ta; - suN::suNgenerator(gen,ta); - std::cout<< "Nc = "< static void testGenerators(void){ - suNmatrix ta; - suNmatrix tb; - std::cout<<"Checking trace ta tb is 0.5 delta_ab"< simd_layout = GridDefaultSimd(4,vComplexF::Nsimd()); - std::vector mpi_layout = GridDefaultMpi(); - std::vector latt_size ({4,4,4,4}); - - GridCartesian Fine(latt_size,simd_layout,mpi_layout); - - LatticeGaugeField Umu(&Fine); + std::vector latt({4,4,4,8}); + GridCartesian * grid = SpaceTimeGrid::makeFourDimGrid(latt, + GridDefaultSimd(Nd,vComplexF::Nsimd()), + GridDefaultMpi()); + + GridRedBlackCartesian * rbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(grid); std::cout<<"*********************************************"<(); - suN::testGenerators<2>(); + SU2::printGenerators(); + SU2::testGenerators(); + std::cout<<"*********************************************"<(); - suN::testGenerators<3>(); - std::cout<<"*********************************************"<(); - suN::testGenerators<4>(); - std::cout<<"*********************************************"<(); - suN::testGenerators<5>(); + SU3::printGenerators(); + SU3::testGenerators(); + + // std::cout<<"*********************************************"< pseeds({1,2,3,4,5}); // once I caught a fish alive + std::vector sseeds({6,7,8,9,10});// then i let it go again + GridParallelRNG pRNG(grid); pRNG.SeedFixedIntegers(pseeds); + GridSerialRNG sRNG; sRNG.SeedFixedIntegers(sseeds); + + // SU3 colour operatoions + LatticeColourMatrix link(grid); + LatticeColourMatrix staple(grid); + int mu=0; + + // Get Staple + ColourWilsonLoops::Staple(staple,Umu,mu); + // Get Link + link = peekIndex(Umu,mu); + + // Apply heatbath to the link + RealD beta=6.0; + int subgroup=0; + int nhb=1; + int trials=0; + int fails=0; + + LatticeInteger one(rbGrid); one = 1; // fill with ones + LatticeInteger mask(grid); mask= zero; + one.checkerboard=Even; + setCheckerboard(mask,one); + + // update Even checkerboard + + SU3::SubGroupHeatBath(sRNG,pRNG,beta,link,staple,subgroup, + nhb,trials,fails,mask); + + Grid_finalize(); }