portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-12 20:27:06 +01:00
Code Releases Activity

Merge branch 'feature/hadrons' of https://github.com/paboyle/Grid into feature/rare_kaon

Browse Source
This commit is contained in:
Lanny91
2017-04-07 11:51:09 +01:00
parent cd1bd921bd ff4e54ef80
commit 5e477ec553
202 changed files with 26985 additions and 2911 deletions
Download Patch File Download Diff File Expand all files Collapse all files

112
tests/IO/Test_nersc_read.cc Normal file
View File

@ -0,0 +1,112 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_nersc_io.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> simd_layout = GridDefaultSimd(4,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
std::vector<int> latt_size = GridDefaultLatt();
int orthodir=3;
int orthosz =latt_size[orthodir];
GridCartesian Fine(latt_size,simd_layout,mpi_layout);
LatticeGaugeField Umu(&Fine);
std::vector<LatticeColourMatrix> U(4,&Fine);
NerscField header;
std::string file("./ckpoint_lat");
NerscIO::readConfiguration(Umu,header,file);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
}
// Painful ; fix syntactical niceness
LatticeComplex LinkTrace(&Fine);
LinkTrace=zero;
for(int mu=0;mu<Nd;mu++){
LinkTrace = LinkTrace + trace(U[mu]);
}
// (1+2+3)=6 = N(N-1)/2 terms
LatticeComplex Plaq(&Fine);
Plaq = zero;
#if 1
for(int mu=1;mu<Nd;mu++){
for(int nu=0;nu<mu;nu++){
Plaq = Plaq + trace(U[mu]*Cshift(U[nu],mu,1)*adj(Cshift(U[mu],nu,1))*adj(U[nu]));
}
}
#endif
double vol = Fine.gSites();
Complex PlaqScale(1.0/vol/6.0/3.0);
std::cout<<GridLogMessage <<"PlaqScale" << PlaqScale<<std::endl;
std::vector<TComplex> Plaq_T(orthosz);
sliceSum(Plaq,Plaq_T,Nd-1);
int Nt = Plaq_T.size();
TComplex Plaq_T_sum;
Plaq_T_sum=zero;
for(int t=0;t<Nt;t++){
Plaq_T_sum = Plaq_T_sum+Plaq_T[t];
Complex Pt=TensorRemove(Plaq_T[t]);
std::cout<<GridLogMessage << "sliced ["<<t<<"]" <<Pt*PlaqScale*Real(Nt)<<std::endl;
}
{
Complex Pt = TensorRemove(Plaq_T_sum);
std::cout<<GridLogMessage << "total " <<Pt*PlaqScale<<std::endl;
}
TComplex Tp = sum(Plaq);
Complex p = TensorRemove(Tp);
std::cout<<GridLogMessage << "calculated plaquettes " <<p*PlaqScale<<std::endl;
Complex LinkTraceScale(1.0/vol/4.0/3.0);
TComplex Tl = sum(LinkTrace);
Complex l = TensorRemove(Tl);
std::cout<<GridLogMessage << "calculated link trace " <<l*LinkTraceScale<<std::endl;
Grid_finalize();
}

3
tests/Test_cayley_even_odd_vec.cc
View File

@ -91,7 +91,8 @@ int main (int argc, char ** argv)
GridParallelRNG sRNG4(sUGrid); sRNG4.SeedFixedIntegers(seeds4);
GridParallelRNG sRNG5(sFGrid); sRNG5.SeedFixedIntegers(seeds5);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4,Umu);
RealD mass=0.1;
RealD M5 =1.8;

18
tests/Test_cshift.cc
View File

@ -61,7 +61,15 @@ int main (int argc, char ** argv)
U=lex;
}
std::stringstream ss;
ss<<"error";
for(int d=0;d<Fine._ndimension;d++){
ss<<"."<<Fine._processor_coor[d];
}
ss<<"_wr_"<<Fine._processor;
std::string fname(ss.str());
std::ofstream ferr(fname);
TComplex cm;
for(int dir=0;dir<4;dir++){
@ -99,11 +107,13 @@ int main (int argc, char ** argv)
Lexicographic::CoorFromIndex(peer,index,latt_size);
if (nrm > 0){
std::cerr<<"FAIL shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cerr<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
ferr<<"FAIL shift "<< shift<<" in dir "<< dir<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
ferr<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
std::cerr<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
ferr<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
}
}}}}
}

211
tests/Test_simd.cc
View File

@ -113,8 +113,6 @@ public:
// outerproduct,
// zeroit
// permute
class funcReduce {
public:
funcReduce() {};
@ -168,7 +166,7 @@ void Tester(const functor &func)
int ok=0;
for(int i=0;i<Nsimd;i++){
if ( abs(reference[i]-result[i])>1.0e-7){
if ( abs(reference[i]-result[i])>1.0e-6){
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< "["<<i<<"] "<< abs(reference[i]-result[i]) << " " <<reference[i]<< " " << result[i]<<std::endl;
ok++;
@ -180,6 +178,65 @@ void Tester(const functor &func)
assert(ok==0);
}
template<class functor>
void IntTester(const functor &func)
{
typedef Integer scal;
typedef vInteger vec;
GridSerialRNG sRNG;
sRNG.SeedRandomDevice();
int Nsimd = vec::Nsimd();
std::vector<scal> input1(Nsimd);
std::vector<scal> input2(Nsimd);
std::vector<scal> result(Nsimd);
std::vector<scal> reference(Nsimd);
std::vector<vec,alignedAllocator<vec> > buf(3);
vec & v_input1 = buf[0];
vec & v_input2 = buf[1];
vec & v_result = buf[2];
for(int i=0;i<Nsimd;i++){
input1[i] = (i + 1) * 30;
input2[i] = (i + 1) * 20;
result[i] = (i + 1) * 10;
}
merge<vec,scal>(v_input1,input1);
merge<vec,scal>(v_input2,input2);
merge<vec,scal>(v_result,result);
func(v_result,v_input1,v_input2);
for(int i=0;i<Nsimd;i++) {
func(reference[i],input1[i],input2[i]);
}
extract<vec,scal>(v_result,result);
std::cout << GridLogMessage << " " << func.name() << std::endl;
std::cout << GridLogDebug << v_input1 << std::endl;
std::cout << GridLogDebug << v_input2 << std::endl;
std::cout << GridLogDebug << v_result << std::endl;
int ok=0;
for(int i=0;i<Nsimd;i++){
if ( reference[i]-result[i] != 0){
std::cout<<GridLogMessage<< "*****" << std::endl;
std::cout<<GridLogMessage<< "["<<i<<"] "<< reference[i]-result[i] << " " <<reference[i]<< " " << result[i]<<std::endl;
ok++;
}
}
if ( ok==0 ) {
std::cout<<GridLogMessage << " OK!" <<std::endl;
}
assert(ok==0);
}
template<class reduced,class scal, class vec,class functor >
void ReductionTester(const functor &func)
@ -245,6 +302,28 @@ public:
}
std::string name(void) const { return std::string("Permute"); }
};
class funcExchange {
public:
int n;
funcExchange(int _n) { n=_n;};
template<class vec> void operator()(vec &r1,vec &r2,vec &i1,vec &i2) const { exchange(r1,r2,i1,i2,n);}
template<class scal> void apply(std::vector<scal> &r1,std::vector<scal> &r2,std::vector<scal> &in1,std::vector<scal> &in2) const {
int sz=in1.size();
int msk = sz>>(n+1);
int j1=0;
int j2=0;
for(int i=0;i<sz;i++) if ( (i&msk) == 0 ) r1[j1++] = in1[ i ];
for(int i=0;i<sz;i++) if ( (i&msk) == 0 ) r1[j1++] = in2[ i ];
for(int i=0;i<sz;i++) if ( (i&msk) ) r2[j2++] = in1[ i ];
for(int i=0;i<sz;i++) if ( (i&msk) ) r2[j2++] = in2[ i ];
}
std::string name(void) const { return std::string("Exchange"); }
};
class funcRotate {
public:
int n;
@ -325,6 +404,89 @@ void PermTester(const functor &func)
assert(ok==0);
}
template<class scal, class vec,class functor >
void ExchangeTester(const functor &func)
{
GridSerialRNG sRNG;
sRNG.SeedRandomDevice();
int Nsimd = vec::Nsimd();
std::vector<scal> input1(Nsimd);
std::vector<scal> input2(Nsimd);
std::vector<scal> result1(Nsimd);
std::vector<scal> result2(Nsimd);
std::vector<scal> reference1(Nsimd);
std::vector<scal> reference2(Nsimd);
std::vector<scal> test1(Nsimd);
std::vector<scal> test2(Nsimd);
std::vector<vec,alignedAllocator<vec> > buf(6);
vec & v_input1 = buf[0];
vec & v_input2 = buf[1];
vec & v_result1 = buf[2];
vec & v_result2 = buf[3];
vec & v_test1 = buf[4];
vec & v_test2 = buf[5];
for(int i=0;i<Nsimd;i++){
random(sRNG,input1[i]);
random(sRNG,input2[i]);
random(sRNG,result1[i]);
random(sRNG,result2[i]);
}
merge<vec,scal>(v_input1,input1);
merge<vec,scal>(v_input2,input2);
merge<vec,scal>(v_result1,result1);
merge<vec,scal>(v_result2,result1);
func(v_result1,v_result2,v_input1,v_input2);
func.apply(reference1,reference2,input1,input2);
func(v_test1,v_test2,v_result1,v_result2);
extract<vec,scal>(v_result1,result1);
extract<vec,scal>(v_result2,result2);
extract<vec,scal>(v_test1,test1);
extract<vec,scal>(v_test2,test2);
std::cout<<GridLogMessage << " " << func.name() << " " <<func.n <<std::endl;
// for(int i=0;i<Nsimd;i++) std::cout << " i "<<i<<" "<<reference1[i]<<" "<<result1[i]<<std::endl;
// for(int i=0;i<Nsimd;i++) std::cout << " i "<<i<<" "<<reference2[i]<<" "<<result2[i]<<std::endl;
for(int i=0;i<Nsimd;i++){
int found=0;
for(int j=0;j<Nsimd;j++){
if(reference1[j]==result1[i]) {
found=1;
// std::cout << " i "<<i<<" j "<<j<<" "<<reference1[j]<<" "<<result1[i]<<std::endl;
}
}
assert(found==1);
}
for(int i=0;i<Nsimd;i++){
int found=0;
for(int j=0;j<Nsimd;j++){
if(reference2[j]==result2[i]) {
found=1;
// std::cout << " i "<<i<<" j "<<j<<" "<<reference2[j]<<" "<<result2[i]<<std::endl;
}
}
assert(found==1);
}
for(int i=0;i<Nsimd;i++){
assert(test1[i]==input1[i]);
assert(test2[i]==input2[i]);
}// std::cout << " i "<< i<<" test1"<<test1[i]<<" "<<input1[i]<<std::endl;
// std::cout << " i "<< i<<" test2"<<test2[i]<<" "<<input2[i]<<std::endl;
// }
}
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
@ -363,6 +525,15 @@ int main (int argc, char ** argv)
PermTester<RealF,vRealF>(funcPermute(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealF exchanges "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
// Log2 iteration
for(int i=0;(1<<i)< vRealF::Nsimd();i++){
ExchangeTester<RealF,vRealF>(funcExchange(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealF rotate "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
@ -394,6 +565,14 @@ int main (int argc, char ** argv)
PermTester<RealD,vRealD>(funcPermute(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealD exchanges "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
// Log2 iteration
for(int i=0;(1<<i)< vRealD::Nsimd();i++){
ExchangeTester<RealD,vRealD>(funcExchange(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vRealD rotate "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
@ -429,6 +608,16 @@ int main (int argc, char ** argv)
PermTester<ComplexF,vComplexF>(funcPermute(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexF exchanges "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
// Log2 iteration
for(int i=0;(1<<i)< vComplexF::Nsimd();i++){
ExchangeTester<ComplexF,vComplexF>(funcExchange(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexF rotate "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
@ -466,12 +655,28 @@ int main (int argc, char ** argv)
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexD exchanges "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
// Log2 iteration
for(int i=0;(1<<i)< vComplexD::Nsimd();i++){
ExchangeTester<ComplexD,vComplexD>(funcExchange(i));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vComplexD rotate "<<std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
for(int r=0;r<vComplexD::Nsimd();r++){
PermTester<ComplexD,vComplexD>(funcRotate(r));
}
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing vInteger "<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
IntTester(funcPlus());
IntTester(funcMinus());
IntTester(funcTimes());
Grid_finalize();
}

23
tests/Test_stencil.cc
View File

@ -66,7 +66,9 @@ int main (int argc, char ** argv)
random(fRNG,Foo);
gaussian(fRNG,Bar);
/*
for (int i=0;i<simd_layout.size();i++){
std::cout <<" simd layout "<<i<<" = "<<simd_layout[i]<<std::endl;
}
Integer stride =1000;
{
double nrm;
@ -79,7 +81,6 @@ int main (int argc, char ** argv)
}
Foo=lex;
}
*/
typedef CartesianStencil<vobj,vobj> Stencil;
for(int dir=0;dir<4;dir++){
@ -92,7 +93,6 @@ int main (int argc, char ** argv)
std::vector<int> displacements(npoint,disp);
Stencil myStencil(&Fine,npoint,0,directions,displacements);
std::vector<int> ocoor(4);
for(int o=0;o<Fine.oSites();o++){
Fine.oCoorFromOindex(ocoor,o);
@ -115,8 +115,11 @@ int main (int argc, char ** argv)
permute(Check._odata[i],Foo._odata[SE->_offset],permute_type);
else if (SE->_is_local)
Check._odata[i] = Foo._odata[SE->_offset];
else
else {
Check._odata[i] = myStencil.CommBuf()[SE->_offset];
// std::cout << " receive "<<i<<" " << Check._odata[i]<<std::endl;
// std::cout << " Foo "<<i<<" " << Foo._odata[i]<<std::endl;
}
}
Real nrmC = norm2(Check);
@ -147,7 +150,12 @@ int main (int argc, char ** argv)
}}}}
if (nrm > 1.0e-4) {
for(int i=0;i<Check._odata.size();i++){
std::cout << i<<" Check.odata "<<Check._odata[i]<< "\n"<<i<<" Bar.odata "<<Bar._odata[i]<<std::endl;
}
}
if (nrm > 1.0e-4) exit(-1);
}
}
@ -182,8 +190,6 @@ int main (int argc, char ** argv)
SimpleCompressor<vobj> compress;
EStencil.HaloExchange(EFoo,compress);
OStencil.HaloExchange(OFoo,compress);
Bar = Cshift(Foo,dir,disp);
@ -196,6 +202,7 @@ int main (int argc, char ** argv)
}
// Implement a stencil code that should agree with that darn cshift!
EStencil.HaloExchange(EFoo,compress);
for(int i=0;i<OCheck._grid->oSites();i++){
int permute_type;
StencilEntry *SE;
@ -209,6 +216,7 @@ int main (int argc, char ** argv)
else
OCheck._odata[i] = EStencil.CommBuf()[SE->_offset];
}
OStencil.HaloExchange(OFoo,compress);
for(int i=0;i<ECheck._grid->oSites();i++){
int permute_type;
StencilEntry *SE;
@ -254,6 +262,7 @@ int main (int argc, char ** argv)
}}}}
if (nrm > 1.0e-4) exit(-1);
}
}

2
tests/core/Test_cf_coarsen_support.cc
View File

@ -64,7 +64,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/core/Test_contfrac_even_odd.cc
View File

@ -70,7 +70,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;

75
tests/core/Test_cshift_red_black.cc
View File

@ -32,6 +32,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
using namespace Grid;
using namespace Grid::QCD;
#define POWER10
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
@ -52,6 +54,7 @@ int main (int argc, char ** argv)
LatticeComplex U(&Fine);
LatticeComplex ShiftU(&Fine);
LatticeComplex rbShiftU(&Fine);
LatticeComplex err(&Fine);
LatticeComplex Ue(&RBFine);
LatticeComplex Uo(&RBFine);
LatticeComplex ShiftUe(&RBFine);
@ -68,7 +71,11 @@ int main (int argc, char ** argv)
Integer i=0;
LatticeCoordinate(coor,d);
lex = lex + coor*stride+i;
#ifndef POWER10
stride=stride*latt_size[d];
#else
stride=stride*10;
#endif
}
U=lex;
}
@ -87,28 +94,31 @@ int main (int argc, char ** argv)
// if ( dir!=1 ) continue;
for(int shift=0;shift<latt_size[dir];shift++){
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir<< " ";
std::cout<<GridLogMessage<<"Even grid"<<std::endl;
// std::cout<<GridLogMessage<<"Even grid"<<std::endl;
ShiftUe = Cshift(Ue,dir,shift); // Shift everything cb by cb
std::cout<<GridLogMessage << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
std::cout<<GridLogMessage<<"Odd grid"<<std::endl;
// std::cout<<GridLogMessage<<"Odd grid"<<std::endl;
ShiftUo = Cshift(Uo,dir,shift);
std::cout<<GridLogMessage << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
std::cout<<GridLogMessage<<"Recombined Even/Odd grids"<<std::endl;
// std::cout<<GridLogMessage<<"Recombined Even/Odd grids"<<std::endl;
setCheckerboard(rbShiftU,ShiftUe);
setCheckerboard(rbShiftU,ShiftUo);
std::cout<<GridLogMessage << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
//std::cout<<GridLogMessage << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
std::cout<<GridLogMessage<<"Full grid shift"<<std::endl;
// std::cout<<GridLogMessage<<"Full grid shift"<<std::endl;
ShiftU = Cshift(U,dir,shift); // Shift everything
std::cout<<GridLogMessage << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
err = ShiftU - rbShiftU;
std::cout<< "\terror " <<norm2(err)<<std::endl;
std::vector<int> coor(4);
std::cout<<GridLogMessage << "Checking the non-checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << " Checking the non-checkerboard shift "<<shift <<" dir "<<dir <<" ... ";
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@ -121,18 +131,26 @@ int main (int argc, char ** argv)
std::vector<int> scoor(coor);
scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
#ifndef POWER10
std::vector<int> powers=latt_size;
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
#else
std::vector<int> powers({1,10,100,1000});
Integer slex = scoor[0]
+ 10 *scoor[1]
+ 100 *scoor[2]
+ 1000 *scoor[3];
#endif
Complex scm(slex);
double nrm = abs(scm-cm()()());
std::vector<int> peer(4);
Complex ctmp = cm;
Integer index=real(ctmp);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
if (nrm > 0){
std::cout<<"FAIL shift "<< shift<<" in dir "<< dir
@ -145,9 +163,10 @@ int main (int argc, char ** argv)
exit(-1);
}
}}}}
std::cout << " OK !"<<std::endl;
int exx=0;
std::cout<<GridLogMessage << "Checking the checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << " Checking the checkerboard shift "<< shift << " dir " << dir <<" ... ";
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@ -157,39 +176,49 @@ int main (int argc, char ** argv)
Integer checkerboard = RBFine.CheckerBoard(coor);
// std::cout << " coor "<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] \n ";
// std::cout << "shift "<< shift <<" dir "<<dir<< " checker board "<< checkerboard << " ";
// std::cout << "Uo " << ShiftUo.checkerboard << " Ue "<<ShiftUe.checkerboard<<std::endl;
if ( checkerboard == ShiftUo.checkerboard ) {
peekSite(cmeo,ShiftUo,coor);
} else {
peekSite(cmeo,ShiftUe,coor);
}
std::vector<int> scoor(coor);
scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
#ifndef POWER10
std::vector<int> powers=latt_size;
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
#else
std::vector<int> powers({1,10,100,1000});
Integer slex = scoor[0]
+ 10 *scoor[1]
+ 100 *scoor[2]
+ 1000 *scoor[3];
#endif
Complex scm(slex);
std::vector<int> peer(4);
Complex ctmp=cmeo;
Integer index=real(ctmp);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
double nrm = abs(cmeo()()()-scm);
if (nrm != 0) {
std::cout << " coor "<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] \n ";
std::cout << "shift "<< shift <<" dir "<<dir<< " checker board "<< checkerboard << " ";
std::cout << "Uo cb = " << ShiftUo.checkerboard << " Ue cb= "<<ShiftUe.checkerboard<<std::endl;
std::cout<<"EOFAIL shift "<< shift<<" in dir "<< dir
<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "
<<" cm " << cm()()()<<" cmeo "
<< cmeo()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
std::cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exx=1;
@ -205,17 +234,17 @@ int main (int argc, char ** argv)
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
std::cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exx=1;
} else if (1) {
} else if (0) {
std::cout<<GridLogMessage<<"PASS shift "<< shift<<" in dir "<< dir
<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
}
}}}}
if (exx) exit(-1);
std::cout << " OK !"<<std::endl;
}
}

74
tests/core/Test_cshift_red_black_rotate.cc
View File

@ -32,6 +32,8 @@ Author: paboyle <paboyle@ph.ed.ac.uk>
using namespace Grid;
using namespace Grid::QCD;
#define POWER10
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
@ -49,6 +51,7 @@ int main (int argc, char ** argv)
GridParallelRNG FineRNG(&Fine); FineRNG.SeedRandomDevice();
LatticeComplex err(&Fine);
LatticeComplex U(&Fine);
LatticeComplex ShiftU(&Fine);
LatticeComplex rbShiftU(&Fine);
@ -66,9 +69,15 @@ int main (int argc, char ** argv)
for(int d=0;d<Nd;d++){
// Integer i=10000;
Integer i=0;
#ifdef POWER10
LatticeCoordinate(coor,Nd-d-1);
lex = lex + coor*stride+i;
stride=stride*10;
#else
LatticeCoordinate(coor,d);
lex = lex + coor*stride+i;
stride=stride*latt_size[d];
#endif
}
U=lex;
}
@ -87,28 +96,30 @@ int main (int argc, char ** argv)
// if ( dir!=1 ) continue;
for(int shift=0;shift<latt_size[dir];shift++){
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir<<std::endl;
std::cout<<GridLogMessage<<"Shifting by "<<shift<<" in direction"<<dir;
std::cout<<GridLogMessage<<"Even grid"<<std::endl;
// std::cout<<GridLogMessage<<"Even grid"<<std::endl;
ShiftUe = Cshift(Ue,dir,shift); // Shift everything cb by cb
std::cout<<GridLogMessage << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
std::cout<<GridLogMessage<<"Odd grid"<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUe " <<norm2(ShiftUe)<<std::endl;
// std::cout<<GridLogMessage<<"Odd grid"<<std::endl;
ShiftUo = Cshift(Uo,dir,shift);
std::cout<<GridLogMessage << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftUo " <<norm2(ShiftUo)<<std::endl;
std::cout<<GridLogMessage<<"Recombined Even/Odd grids"<<std::endl;
// std::cout<<GridLogMessage<<"Recombined Even/Odd grids"<<std::endl;
setCheckerboard(rbShiftU,ShiftUe);
setCheckerboard(rbShiftU,ShiftUo);
std::cout<<GridLogMessage << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<GridLogMessage << "\trbShiftU " <<norm2(rbShiftU)<<std::endl;
std::cout<<GridLogMessage<<"Full grid shift"<<std::endl;
// std::cout<<GridLogMessage<<"Full grid shift"<<std::endl;
ShiftU = Cshift(U,dir,shift); // Shift everything
std::cout<<GridLogMessage << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
// std::cout<<GridLogMessage << "\tShiftU " <<norm2(rbShiftU)<<std::endl;
err = ShiftU - rbShiftU;
std::cout<< "\terror " <<norm2(err)<<std::endl;
std::vector<int> coor(4);
std::cout<<GridLogMessage << "Checking the non-checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << " Checking the non-checkerboard shift "<< shift << " dir "<<dir <<"... ";
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@ -120,11 +131,20 @@ int main (int argc, char ** argv)
std::vector<int> scoor(coor);
scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
#ifdef POWER10
std::vector<int> powers({1,10,100,1000});
Integer slex = scoor[3]
+ 10 *scoor[2]
+ 100 *scoor[1]
+ 1000 *scoor[0];
#else
std::vector<int> powers=latt_size;
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
#endif
Complex scm(slex);
@ -132,7 +152,7 @@ int main (int argc, char ** argv)
std::vector<int> peer(4);
Complex ctmp = cm;
Integer index=real(ctmp);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
if (nrm > 0){
std::cout<<"FAIL shift "<< shift<<" in dir "<< dir
@ -140,14 +160,16 @@ int main (int argc, char ** argv)
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
std::cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exit(-1);
}
}}}}
std::cout<< " OK! "<<std::endl ;
int exx=0;
std::cout<<GridLogMessage << "Checking the checkerboard shift"<<std::endl;
std::cout<<GridLogMessage << " Checking the checkerboard shift "<< shift << " dir " << dir <<"... ";
for(coor[3]=0;coor[3]<latt_size[3];coor[3]++){
for(coor[2]=0;coor[2]<latt_size[2];coor[2]++){
for(coor[1]=0;coor[1]<latt_size[1];coor[1]++){
@ -169,18 +191,26 @@ int main (int argc, char ** argv)
std::vector<int> scoor(coor);
scoor[dir] = (scoor[dir]+shift)%latt_size[dir];
#ifdef POWER10
std::vector<int> powers({1,10,100,1000});
Integer slex = scoor[3]
+ 10 *scoor[2]
+ 100 *scoor[1]
+ 1000 *scoor[0];
#else
std::vector<int> powers = latt_size;
Integer slex = scoor[0]
+ latt_size[0]*scoor[1]
+ latt_size[0]*latt_size[1]*scoor[2]
+ latt_size[0]*latt_size[1]*latt_size[2]*scoor[3];
#endif
Complex scm(slex);
std::vector<int> peer(4);
Complex ctmp=cmeo;
Integer index=real(ctmp);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
double nrm = abs(cmeo()()()-scm);
if (nrm != 0) {
@ -189,10 +219,9 @@ int main (int argc, char ** argv)
<< cmeo()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
std::cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exx=1;
}
ctmp=cm;
@ -205,16 +234,17 @@ int main (int argc, char ** argv)
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
std::cout<<"Got "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
index=real(scm);
Lexicographic::CoorFromIndex(peer,index,latt_size);
Lexicographic::CoorFromIndex(peer,index,powers);
std::cout<<"Expect "<<index<<" " << peer[0]<<","<<peer[1]<<","<<peer[2]<<","<<peer[3]<<std::endl;
exx=1;
} else if (1) {
} else if (0) {
std::cout<<GridLogMessage<<"PASS shift "<< shift<<" in dir "<< dir
<<" ["<<coor[0]<<","<<coor[1]<<","<<coor[2]<<","<<coor[3]<<"] = "
<< cm()()()<<" expect "<<scm<<" "<<nrm<<std::endl;
}
}}}}
if (exx) exit(-1);
std::cout<< " OK! "<<std::endl ;
}
}

1
tests/core/Test_cshift_rotate.cc
View File

@ -45,7 +45,6 @@ int main (int argc, char ** argv)
LatticeComplex U(&Fine);
LatticeComplex ShiftU(&Fine);
LatticeComplex lex(&Fine);
lex=zero;
Integer stride =1;

2
tests/core/Test_dwf_even_odd.cc
View File

@ -72,7 +72,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid); tmp=zero;
LatticeFermion err(FGrid); tmp=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)

2
tests/core/Test_dwf_rb5d.cc
View File

@ -81,7 +81,7 @@ int main (int argc, char ** argv)
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)

2
tests/core/Test_gpwilson_even_odd.cc
View File

@ -61,7 +61,7 @@ int main (int argc, char ** argv)
FermionField ref(&Grid); ref=zero;
FermionField tmp(&Grid); tmp=zero;
FermionField err(&Grid); tmp=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

291
tests/core/Test_staggered.cc Normal file
View File

@ -0,0 +1,291 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_wilson.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
GridCartesian Grid(latt_size,simd_layout,mpi_layout);
GridRedBlackCartesian RBGrid(latt_size,simd_layout,mpi_layout);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(seeds);
// pRNG.SeedRandomDevice();
typedef typename ImprovedStaggeredFermionR::FermionField FermionField;
typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField;
typename ImprovedStaggeredFermionR::ImplParams params;
FermionField src (&Grid); random(pRNG,src);
FermionField result(&Grid); result=zero;
FermionField ref(&Grid); ref=zero;
FermionField tmp(&Grid); tmp=zero;
FermionField err(&Grid); tmp=zero;
FermionField phi (&Grid); random(pRNG,phi);
FermionField chi (&Grid); random(pRNG,chi);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;
for(int mu=0;mu<Nd;mu++){
volume=volume*latt_size[mu];
}
// Only one non-zero (y)
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu,mu);
/* Debug force unit
U[mu] = 1.0;
PokeIndex<LorentzIndex>(Umu,U[mu],mu);
*/
}
ref = zero;
RealD mass=0.1;
RealD c1=9.0/8.0;
RealD c2=-1.0/24.0;
RealD u0=1.0;
{ // Simple improved staggered implementation
ref = zero;
RealD c1tad = 0.5*c1/u0;
RealD c2tad = 0.5*c2/u0/u0/u0;
Lattice<iScalar<vInteger> > coor(&Grid);
Lattice<iScalar<vInteger> > x(&Grid); LatticeCoordinate(x,0);
Lattice<iScalar<vInteger> > y(&Grid); LatticeCoordinate(y,1);
Lattice<iScalar<vInteger> > z(&Grid); LatticeCoordinate(z,2);
Lattice<iScalar<vInteger> > t(&Grid); LatticeCoordinate(t,3);
Lattice<iScalar<vInteger> > lin_z(&Grid); lin_z=x+y;
Lattice<iScalar<vInteger> > lin_t(&Grid); lin_t=x+y+z;
for(int mu=0;mu<Nd;mu++){
// Staggered Phase.
ComplexField phases(&Grid); phases=1.0;
if ( mu == 1 ) phases = where( mod(x ,2)==(Integer)0, phases,-phases);
if ( mu == 2 ) phases = where( mod(lin_z,2)==(Integer)0, phases,-phases);
if ( mu == 3 ) phases = where( mod(lin_t,2)==(Integer)0, phases,-phases);
tmp = PeriodicBC::CovShiftForward(U[mu],mu,src);
ref = ref +c1tad*tmp*phases; // Forward 1 hop
tmp = PeriodicBC::CovShiftForward(U[mu],mu,tmp); // 2 hop
tmp = PeriodicBC::CovShiftForward(U[mu],mu,tmp); // 3 hop
ref = ref +c2tad*tmp*phases; // Forward 3 hop
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,src);
ref = ref -c1tad*tmp*phases; // Forward 3 hop
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,tmp);
tmp = PeriodicBC::CovShiftBackward(U[mu],mu,tmp);
ref = ref -c2tad*tmp*phases; // Forward 3 hop
}
// ref = ref + mass * src;
}
ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass,c1,c2,u0,params);
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
int ncall=1000;
double t0=usecond();
for(int i=0;i<ncall;i++){
Ds.Dhop(src,result,0);
}
double t1=usecond();
double t2;
double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + == 1146
std::cout<<GridLogMessage << "Called Ds"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
err = ref-result;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
FermionField src_e (&RBGrid);
FermionField src_o (&RBGrid);
FermionField r_e (&RBGrid);
FermionField r_o (&RBGrid);
FermionField r_eo (&Grid);
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ds.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Ds.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
Ds.Dhop (src,ref,DaggerNo);
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
err= ref - r_eo;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
FermionField chi_e (&RBGrid);
FermionField chi_o (&RBGrid);
FermionField dchi_e (&RBGrid);
FermionField dchi_o (&RBGrid);
FermionField phi_e (&RBGrid);
FermionField phi_o (&RBGrid);
FermionField dphi_e (&RBGrid);
FermionField dphi_o (&RBGrid);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
Ds.Meooe(chi_e,dchi_o);
Ds.Meooe(chi_o,dchi_e);
Ds.MeooeDag(phi_e,dphi_o);
Ds.MeooeDag(phi_o,dphi_e);
ComplexD pDce = innerProduct(phi_e,dchi_e);
ComplexD pDco = innerProduct(phi_o,dchi_o);
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ds.Mooee(chi_e,src_e);
Ds.MooeeInv(src_e,phi_e);
Ds.Mooee(chi_o,src_o);
Ds.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ds.MooeeDag(chi_e,src_e);
Ds.MooeeInvDag(src_e,phi_e);
Ds.MooeeDag(chi_o,src_o);
Ds.MooeeInvDag(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(pRNG,phi);
random(pRNG,chi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<ImprovedStaggeredFermionR,FermionField> HermOpEO(Ds);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}

309
tests/core/Test_staggered5D.cc Normal file
View File

@ -0,0 +1,309 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_wilson.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
std::cout << GridLogMessage << "Making s innermost grids"<<std::endl;
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALF"<< sizeof(RealF)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REALD"<< sizeof(RealD)<<std::endl;
std::cout<<GridLogMessage << "Grid floating point word size is REAL"<< sizeof(Real)<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG4(UGrid);
GridParallelRNG pRNG5(FGrid);
pRNG4.SeedFixedIntegers(seeds);
pRNG5.SeedFixedIntegers(seeds);
typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField;
typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField;
typename ImprovedStaggeredFermion5DR::ImplParams params;
FermionField src (FGrid);
random(pRNG5,src);
FermionField result(FGrid); result=zero;
FermionField ref(FGrid); ref=zero;
FermionField tmp(FGrid); tmp=zero;
FermionField err(FGrid); tmp=zero;
FermionField phi (FGrid); random(pRNG5,phi);
FermionField chi (FGrid); random(pRNG5,chi);
LatticeGaugeField Umu(UGrid); SU3::ColdConfiguration(pRNG4,Umu);
LatticeGaugeField Umua(UGrid); Umua=Umu;
double volume=Ls;
for(int mu=0;mu<Nd;mu++){
volume=volume*latt_size[mu];
}
////////////////////////////////////
// Naive implementation needs to
// replicate across fifth dimension
////////////////////////////////////
LatticeGaugeField Umu5d(FGrid);
for(int ss=0;ss<Umu._grid->oSites();ss++){
for(int s=0;s<Ls;s++){
Umu5d._odata[Ls*ss+s] = Umu._odata[ss];
}
}
std::vector<LatticeColourMatrix> U(4,FGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = PeekIndex<LorentzIndex>(Umu5d,mu);
}
RealD mass=0.1;
RealD c1=9.0/8.0;
RealD c2=-1.0/24.0;
RealD u0=1.0;
{ // Simple improved staggered implementation
ref = zero;
RealD c1tad = 0.5*c1/u0;
RealD c2tad = 0.5*c2/u0/u0/u0;
Lattice<iScalar<vInteger> > coor(FGrid);
Lattice<iScalar<vInteger> > x(FGrid); LatticeCoordinate(x,1); // s innermost
Lattice<iScalar<vInteger> > y(FGrid); LatticeCoordinate(y,2);
Lattice<iScalar<vInteger> > z(FGrid); LatticeCoordinate(z,3);
Lattice<iScalar<vInteger> > t(FGrid); LatticeCoordinate(t,4);
Lattice<iScalar<vInteger> > lin_z(FGrid); lin_z=x+y;
Lattice<iScalar<vInteger> > lin_t(FGrid); lin_t=x+y+z;
for(int mu=0;mu<Nd;mu++){
// Staggered Phase.
ComplexField phases(FGrid); phases=1.0;
if ( mu == 1 ) phases = where( mod(x ,2)==(Integer)0, phases,-phases);
if ( mu == 2 ) phases = where( mod(lin_z,2)==(Integer)0, phases,-phases);
if ( mu == 3 ) phases = where( mod(lin_t,2)==(Integer)0, phases,-phases);
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,src);
ref = ref +c1tad*tmp*phases; // Forward 1 hop
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,tmp); // 2 hop
tmp = PeriodicBC::CovShiftForward(U[mu],mu+1,tmp); // 3 hop
ref = ref +c2tad*tmp*phases; // Forward 3 hop
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,src);
ref = ref -c1tad*tmp*phases; // Forward 3 hop
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,tmp);
tmp = PeriodicBC::CovShiftBackward(U[mu],mu+1,tmp);
ref = ref -c2tad*tmp*phases; // Forward 3 hop
}
}
ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,c1,c2,u0,params);
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage << "Calling Ds"<<std::endl;
int ncall=1000;
double t0=usecond();
for(int i=0;i<ncall;i++){
Ds.Dhop(src,result,0);
}
double t1=usecond();
double t2;
double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + == 1146
std::cout<<GridLogMessage << "Called Ds"<<std::endl;
// std::cout<<GridLogMessage << "result"<< result <<std::endl;
// std::cout<<GridLogMessage << "ref "<< ref <<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "norm ref "<< norm2(ref)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
err = ref-result;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Deo + Doe = Dunprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
FermionField src_e (FrbGrid);
FermionField src_o (FrbGrid);
FermionField r_e (FrbGrid);
FermionField r_o (FrbGrid);
FermionField r_eo (FGrid);
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ds.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Ds.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
Ds.Dhop (src,ref,DaggerNo);
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
err= ref - r_eo;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
FermionField chi_e (FrbGrid);
FermionField chi_o (FrbGrid);
FermionField dchi_e (FrbGrid);
FermionField dchi_o (FrbGrid);
FermionField phi_e (FrbGrid);
FermionField phi_o (FrbGrid);
FermionField dphi_e (FrbGrid);
FermionField dphi_o (FrbGrid);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
Ds.Meooe(chi_e,dchi_o);
Ds.Meooe(chi_o,dchi_e);
Ds.MeooeDag(phi_e,dphi_o);
Ds.MeooeDag(phi_o,dphi_e);
ComplexD pDce = innerProduct(phi_e,dchi_e);
ComplexD pDco = innerProduct(phi_o,dchi_o);
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ds.Mooee(chi_e,src_e);
Ds.MooeeInv(src_e,phi_e);
Ds.Mooee(chi_o,src_o);
Ds.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ds.MooeeDag(chi_e,src_e);
Ds.MooeeInvDag(src_e,phi_e);
Ds.MooeeDag(chi_o,src_o);
Ds.MooeeInvDag(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(pRNG5,phi);
random(pRNG5,chi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
SchurDiagMooeeOperator<ImprovedStaggeredFermion5DR,FermionField> HermOpEO(Ds);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}

193
tests/core/Test_staggered5Dvec.cc Normal file
View File

@ -0,0 +1,193 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./benchmarks/Benchmark_wilson.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
std::vector<int> latt_size = GridDefaultLatt();
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplex::Nsimd());
std::vector<int> mpi_layout = GridDefaultMpi();
const int Ls=16;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::cout << GridLogMessage << "Making s innermost grids"<<std::endl;
GridCartesian * sUGrid = SpaceTimeGrid::makeFourDimDWFGrid(GridDefaultLatt(),GridDefaultMpi());
GridRedBlackCartesian * sUrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(sUGrid);
GridCartesian * sFGrid = SpaceTimeGrid::makeFiveDimDWFGrid(Ls,UGrid);
GridRedBlackCartesian * sFrbGrid = SpaceTimeGrid::makeFiveDimDWFRedBlackGrid(Ls,UGrid);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG4(UGrid);
GridParallelRNG pRNG5(FGrid);
pRNG4.SeedFixedIntegers(seeds);
pRNG5.SeedFixedIntegers(seeds);
typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField;
typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField;
typename ImprovedStaggeredFermion5DR::ImplParams params;
FermionField src (FGrid);
random(pRNG5,src);
/*
std::vector<int> site({0,1,2,0,0});
ColourVector cv = zero;
cv()()(0)=1.0;
src = zero;
pokeSite(cv,src,site);
*/
FermionField result(FGrid); result=zero;
FermionField tmp(FGrid); tmp=zero;
FermionField err(FGrid); tmp=zero;
FermionField phi (FGrid); random(pRNG5,phi);
FermionField chi (FGrid); random(pRNG5,chi);
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(pRNG4,Umu);
/*
for(int mu=1;mu<4;mu++){
auto tmp = PeekIndex<LorentzIndex>(Umu,mu);
tmp = zero;
PokeIndex<LorentzIndex>(Umu,tmp,mu);
}
*/
double volume=Ls;
for(int mu=0;mu<Nd;mu++){
volume=volume*latt_size[mu];
}
RealD mass=0.1;
RealD c1=9.0/8.0;
RealD c2=-1.0/24.0;
RealD u0=1.0;
ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,c1,c2,u0,params);
ImprovedStaggeredFermionVec5dR sDs(Umu,Umu,*sFGrid,*sFrbGrid,*sUGrid,*sUrbGrid,mass,c1,c2,u0,params);
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing Dhop against cshift implementation "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
int ncall=1000;
int ncall1=1000;
double t0(0),t1(0);
double flops=(16*(3*(6+8+8)) + 15*3*2)*volume*ncall; // == 66*16 + == 1146
std::cout<<GridLogMessage << "Calling staggered operator"<<std::endl;
t0=usecond();
for(int i=0;i<ncall1;i++){
Ds.Dhop(src,result,0);
}
t1=usecond();
std::cout<<GridLogMessage << "Called Ds"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(result)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
std::cout<<GridLogMessage << "Calling vectorised staggered operator"<<std::endl;
#ifdef AVX512
QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptInlineAsm;
#else
QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptGeneric;
#endif
t0=usecond();
for(int i=0;i<ncall1;i++){
Ds.Dhop(src,tmp,0);
}
t1=usecond();
std::cout<<GridLogMessage << "Called Ds ASM"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(tmp)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
err = tmp-result;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
FermionField ssrc (sFGrid); localConvert(src,ssrc);
FermionField sresult(sFGrid); sresult=zero;
QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptHandUnroll;
t0=usecond();
for(int i=0;i<ncall1;i++){
sDs.Dhop(ssrc,sresult,0);
}
t1=usecond();
localConvert(sresult,tmp);
std::cout<<GridLogMessage << "Called sDs unroll"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(sresult)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)<<std::endl;
#ifdef AVX512
QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptInlineAsm;
#else
QCD::StaggeredKernelsStatic::Opt=QCD::StaggeredKernelsStatic::OptGeneric;
#endif
err = tmp-result;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
int extra=1;
t0=usecond();
for(int i=0;i<ncall1*extra;i++){
sDs.Dhop(ssrc,sresult,0);
}
t1=usecond();
localConvert(sresult,tmp);
std::cout<<GridLogMessage << "Called sDs asm"<<std::endl;
std::cout<<GridLogMessage << "norm result "<< norm2(sresult)<<std::endl;
std::cout<<GridLogMessage << "mflop/s = "<< flops/(t1-t0)*extra<<std::endl;
err = tmp-result;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<<std::endl;
Grid_finalize();
}

2
tests/core/Test_wilson_even_odd.cc
View File

@ -71,7 +71,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(&Grid); ref=zero;
LatticeFermion tmp(&Grid); tmp=zero;
LatticeFermion err(&Grid); tmp=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

2
tests/core/Test_wilson_tm_even_odd.cc
View File

@ -70,7 +70,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(&Grid); ref=zero;
LatticeFermion tmp(&Grid); tmp=zero;
LatticeFermion err(&Grid); tmp=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);
double volume=1;

273
tests/core/Test_zmobius_even_odd.cc Normal file
View File

@ -0,0 +1,273 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_even_odd.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template<class d>
struct scal {
d internal;
};
Gamma::Algebra Gmu [] = {
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT
};
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
const int Ls=8;
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplex::Nsimd()),GridDefaultMpi());
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
std::vector<int> seeds4({1,2,3,4});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
LatticeFermion src (FGrid); random(RNG5,src);
LatticeFermion phi (FGrid); random(RNG5,phi);
LatticeFermion chi (FGrid); random(RNG5,chi);
LatticeFermion result(FGrid); result=zero;
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid); tmp=zero;
LatticeFermion err(FGrid); tmp=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
// Only one non-zero (y)
Umu=zero;
for(int nn=0;nn<Nd;nn++){
random(RNG4,U[nn]);
if ( nn>0 )
U[nn]=zero;
PokeIndex<LorentzIndex>(Umu,U[nn],nn);
}
RealD mass=0.1;
RealD M5 =1.8;
std::vector < std::complex<double> > omegas;
for(int i=0;i<Ls;i++){
double imag = 0.;
if (i==0) imag=1.;
if (i==Ls-1) imag=-1.;
std::complex<double> temp (0.25+0.01*i, imag*0.1);
omegas.push_back(temp);
}
ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,1.,0.);
// DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
LatticeFermion src_e (FrbGrid);
LatticeFermion src_o (FrbGrid);
LatticeFermion r_e (FrbGrid);
LatticeFermion r_o (FrbGrid);
LatticeFermion r_eo (FGrid);
LatticeFermion r_eeoo(FGrid);
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Testing that Meo + Moe + Moo + Mee = Munprec "<<std::endl;
std::cout<<GridLogMessage<<"=========================================================="<<std::endl;
pickCheckerboard(Even,src_e,src);
pickCheckerboard(Odd,src_o,src);
Ddwf.Meooe(src_e,r_o); std::cout<<GridLogMessage<<"Applied Meo"<<std::endl;
Ddwf.Meooe(src_o,r_e); std::cout<<GridLogMessage<<"Applied Moe"<<std::endl;
setCheckerboard(r_eo,r_o);
setCheckerboard(r_eo,r_e);
Ddwf.Mooee(src_e,r_e); std::cout<<GridLogMessage<<"Applied Mee"<<std::endl;
Ddwf.Mooee(src_o,r_o); std::cout<<GridLogMessage<<"Applied Moo"<<std::endl;
setCheckerboard(r_eeoo,r_e);
setCheckerboard(r_eeoo,r_o);
r_eo=r_eo+r_eeoo;
Ddwf.M(src,ref);
// std::cout<<GridLogMessage << r_eo<<std::endl;
// std::cout<<GridLogMessage << ref <<std::endl;
err= ref - r_eo;
std::cout<<GridLogMessage << "EO norm diff "<< norm2(err)<< " "<<norm2(ref)<< " " << norm2(r_eo) <<std::endl;
LatticeComplex cerr(FGrid);
cerr = localInnerProduct(err,err);
// std::cout<<GridLogMessage << cerr<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MooeeDagger is the dagger of Mooee by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
LatticeFermion chi_e (FrbGrid);
LatticeFermion chi_o (FrbGrid);
LatticeFermion dchi_e (FrbGrid);
LatticeFermion dchi_o (FrbGrid);
LatticeFermion phi_e (FrbGrid);
LatticeFermion phi_o (FrbGrid);
LatticeFermion dphi_e (FrbGrid);
LatticeFermion dphi_o (FrbGrid);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
Ddwf.Mooee(chi_e,dchi_o);
Ddwf.Mooee(chi_o,dchi_e);
Ddwf.MooeeDag(phi_e,dphi_o);
Ddwf.MooeeDag(phi_o,dphi_e);
ComplexD pDce = innerProduct(phi_e,dchi_e);
ComplexD pDco = innerProduct(phi_o,dchi_o);
ComplexD cDpe = innerProduct(chi_e,dphi_e);
ComplexD cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test Ddagger is the dagger of D by requiring "<<std::endl;
std::cout<<GridLogMessage<<"= < phi | Deo | chi > * = < chi | Deo^dag| phi> "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
Ddwf.Meooe(chi_e,dchi_o);
Ddwf.Meooe(chi_o,dchi_e);
Ddwf.MeooeDag(phi_e,dphi_o);
Ddwf.MeooeDag(phi_o,dphi_e);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDce-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDco-conj(cDpe) <<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInv Mee = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ddwf.Mooee(chi_e,src_e);
Ddwf.MooeeInv(src_e,phi_e);
Ddwf.Mooee(chi_o,src_o);
Ddwf.MooeeInv(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MeeInvDag MeeDag = 1 "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
Ddwf.MooeeDag(chi_e,src_e);
Ddwf.MooeeInvDag(src_e,phi_e);
Ddwf.MooeeDag(chi_o,src_o);
Ddwf.MooeeInvDag(src_o,phi_o);
setCheckerboard(phi,phi_e);
setCheckerboard(phi,phi_o);
err = phi-chi;
std::cout<<GridLogMessage << "norm diff "<< norm2(err)<< std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
std::cout<<GridLogMessage<<"= Test MpcDagMpc is Hermitian "<<std::endl;
std::cout<<GridLogMessage<<"=============================================================="<<std::endl;
random(RNG5,phi);
random(RNG5,chi);
pickCheckerboard(Even,chi_e,chi);
pickCheckerboard(Odd ,chi_o,chi);
pickCheckerboard(Even,phi_e,phi);
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
SchurDiagMooeeOperator<ZMobiusFermionR,LatticeFermion> HermOpEO(Ddwf);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);
HermOpEO.MpcDagMpc(phi_e,dphi_e,t1,t2);
HermOpEO.MpcDagMpc(phi_o,dphi_o,t1,t2);
pDce = innerProduct(phi_e,dchi_e);
pDco = innerProduct(phi_o,dchi_o);
cDpe = innerProduct(chi_e,dphi_e);
cDpo = innerProduct(chi_o,dphi_o);
std::cout<<GridLogMessage <<"e "<<pDce<<" "<<cDpe <<std::endl;
std::cout<<GridLogMessage <<"o "<<pDco<<" "<<cDpo <<std::endl;
std::cout<<GridLogMessage <<"pDce - conj(cDpo) "<< pDco-conj(cDpo) <<std::endl;
std::cout<<GridLogMessage <<"pDco - conj(cDpe) "<< pDce-conj(cDpe) <<std::endl;
Grid_finalize();
}

6
tests/debug/Test_cayley_cg.cc
View File

@ -106,10 +106,16 @@ int main (int argc, char ** argv)
RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5;
std::vector<ComplexD> gamma(Ls,ComplexD(1.0,0.0));
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestCGinversions<MobiusFermionR>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c);
TestCGinversions<ZMobiusFermionR>(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
TestCGinversions<MobiusZolotarevFermionR>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);

2
tests/debug/Test_cayley_coarsen_support.cc
View File

@ -77,7 +77,7 @@ int main (int argc, char ** argv)
LatticeFermion ref(FGrid); ref=zero;
LatticeFermion tmp(FGrid);
LatticeFermion err(FGrid);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
#if 0
std::vector<LatticeColourMatrix> U(4,UGrid);

8
tests/debug/Test_cayley_even_odd.cc
View File

@ -70,7 +70,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;
@ -81,10 +81,16 @@ int main (int argc, char ** argv)
RealD b=1.5;// Scale factor b+c=2, b-c=1
RealD c=0.5;
std::vector<ComplexD> gamma(Ls,ComplexD(1.0,0.1));
std::cout<<GridLogMessage <<"MobiusFermion test"<<std::endl;
MobiusFermionR Dmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
TestWhat<MobiusFermionR>(Dmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"ZMobiusFermion test"<<std::endl;
ZMobiusFermionR ZDmob(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,gamma,b,c);
TestWhat<ZMobiusFermionR>(ZDmob,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);
std::cout<<GridLogMessage <<"MobiusZolotarevFermion test"<<std::endl;
MobiusZolotarevFermionR Dzolo(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c,0.1,2.0);
TestWhat<MobiusZolotarevFermionR>(Dzolo,FGrid,FrbGrid,UGrid,UrbGrid,mass,M5,&RNG4,&RNG5);

3
tests/debug/Test_zmm.cc
View File

@ -26,7 +26,6 @@ See the full license in the file "LICENSE" in the top level distribution directo
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
#include <Grid/PerfCount.h>
#ifdef TEST_ZMM
@ -187,7 +186,7 @@ int main(int argc,char **argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
random(RNG5,src);
#if 1
random(RNG4,Umu);
SU3::HotConfiguration(RNG4,Umu);
#else
int mmu=2;
std::vector<LatticeColourMatrix> U(4,UGrid);

2
tests/forces/Test_contfrac_force.cc
View File

@ -31,8 +31,6 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);

2
tests/forces/Test_dwf_force.cc
View File

@ -31,8 +31,6 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{
Grid_init(&argc,&argv);

2
tests/forces/Test_dwf_gpforce.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_gp_rect_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_gpdwf_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_gpwilson_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_partfrac_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_rect_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_wilson_force.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_wilson_force_phiMdagMphi.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

2
tests/forces/Test_wilson_force_phiMphi.cc
View File

@ -31,7 +31,7 @@ using namespace std;
using namespace Grid;
using namespace Grid::QCD;
#define parallel_for PARALLEL_FOR_LOOP for
int main (int argc, char ** argv)
{

78
tests/qdpxx/Test_qdpxx_loops_staples.cc
View File

@ -282,8 +282,8 @@ double calc_grid_p(Grid::QCD::LatticeGaugeField & Umu)
Grid::QCD::LatticeColourMatrix tmp(UGrid);
tmp = Grid::zero;
Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,2);
Grid::QCD::PokeIndex<Grid::QCD::LorentzIndex>(Umu,tmp,3);
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,2);
Grid::QCD::PokeIndex<LorentzIndex>(Umu,tmp,3);
Grid::QCD::WilsonGaugeActionR Wilson(beta); // Just take beta = 1.0
@ -311,7 +311,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
std::vector<Grid::QCD::LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){
U[mu] = Grid::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
U[mu] = Grid::PeekIndex<LorentzIndex>(Umu,mu);
}
Grid::QCD::LatticeComplex rect(UGrid);
@ -322,7 +322,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
for(int nu=0;nu<Grid::QCD::Nd;nu++){
if ( mu!=nu ) {
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::traceDirRectangle(rect,U,mu,nu);
Grid::QCD::ColourWilsonLoops::traceDirRectangle(rect,U,mu,nu);
trect = Grid::sum(rect);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
@ -344,10 +344,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ ___
// | __ |
Stap =
Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
Grid::QCD::CovShiftForward (U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
@ -361,10 +361,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __
// |__ __ |
Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu, U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
@ -375,10 +375,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __
// |__ __ |
Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
@ -390,10 +390,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// __ ___
// |__ |
Stap = Grid::Cshift(Grid::QCD::CovShiftForward (U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,U[mu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
@ -412,12 +412,12 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
* Make staple for loops centered at coor of link ; this one is ok. // |
*/
// Stap =
// Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
// Grid::adj(Grid::QCD::CovShiftForward(U[nu],nu,Grid::QCD::CovShiftForward(U[nu],nu,U[mu]))) ;
Stap = Grid::Cshift(Grid::QCD::CovShiftForward(U[nu],nu,
Grid::QCD::CovShiftForward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
// Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]),mu,1)* // ->||
// Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[mu]))) ;
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu, Grid::Cshift(adj(U[nu]),nu,-1))))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
SumTrStap += TrStap;
@ -433,10 +433,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// | |
// --
Stap = Grid::Cshift(Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[nu],nu,
Grid::QCD::CovShiftBackward(U[mu],mu,
Grid::QCD::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
Stap = Grid::Cshift(Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[nu],nu,
Grid::QCD::PeriodicBC::CovShiftBackward(U[mu],mu,
Grid::QCD::PeriodicBC::CovShiftForward (U[nu],nu,U[nu])))) , mu, 1);
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
@ -460,10 +460,10 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
Grid::QCD::LatticeColourMatrix tmp(UGrid);
// 2 (mu)x1(nu)
left_2= Grid::QCD::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
left_2= Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,U[mu]); // Umu(x) Umu(x+mu)
tmp=Grid::Cshift(U[nu],mu,2); // Unu(x+2mu)
upper_l= Grid::QCD::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(tmp,nu,Grid::adj(left_2)); // Unu(x+2mu) Umu^dag(x+mu+nu) Umu^dag(x+nu)
// __ __
// = |
@ -533,9 +533,9 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
// _
// | |
// | |
Grid::QCD::LatticeColourMatrix up2= Grid::QCD::CovShiftForward(U[nu],nu,U[nu]);
Grid::QCD::LatticeColourMatrix up2= Grid::QCD::PeriodicBC::CovShiftForward(U[nu],nu,U[nu]);
upper_l= Grid::QCD::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
upper_l= Grid::QCD::PeriodicBC::CovShiftForward(Grid::Cshift(up2,mu,1),nu,Grid::Cshift(adj(U[mu]),nu,1));
ds_U= upper_l*Grid::adj(up2);
RectPlaq_d = Grid::trace(U[mu]*ds_U);
@ -555,7 +555,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
downer_l= |
(x)<----V
*/
down_l= Grid::adj(Grid::QCD::CovShiftForward(U[mu],mu,up2)); //downer_l
down_l= Grid::adj(Grid::QCD::PeriodicBC::CovShiftForward(U[mu],mu,up2)); //downer_l
/*
^ |
down_staple = | V
@ -616,9 +616,9 @@ void check_grid_r_staple(Grid::QCD::LatticeGaugeField & Umu)
// Vol as for each site
Grid::RealD RectScale(1.0/vol/12.0/6.0/3.0);
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::RectStaple(staple,Umu,mu);
Grid::QCD::ColourWilsonLoops::RectStaple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * RectScale;
Grid::QCD::TComplex Tp = Grid::sum(Traced);
@ -655,9 +655,9 @@ void check_grid_p_staple(Grid::QCD::LatticeGaugeField & Umu)
// Vol as for each site
Grid::RealD Scale(1.0/vol/12.0/2.0/3.0);
Grid::QCD::WilsonLoops<Grid::QCD::LatticeGaugeField>::Staple(staple,Umu,mu);
Grid::QCD::ColourWilsonLoops::Staple(staple,Umu,mu);
link = Grid::QCD::PeekIndex<Grid::QCD::LorentzIndex>(Umu,mu);
link = Grid::QCD::PeekIndex<LorentzIndex>(Umu,mu);
Traced = Grid::trace( link*staple) * Scale;
Grid::QCD::TComplex Tp = Grid::sum(Traced);

364
tests/qdpxx/Test_qdpxx_stag.cc Normal file
View File

@ -0,0 +1,364 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/qdpxx/Test_qdpxx_munprec.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
Author: paboyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
double mq=0.1;
typedef Grid::QCD::StaggeredImplR::FermionField FermionField;
typedef Grid::QCD::LatticeGaugeField GaugeField;
void make_gauge (GaugeField & lat, FermionField &src);
void calc_grid (GaugeField & lat, GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
void calc_chroma (GaugeField & lat,GaugeField & uthin,GaugeField & ufat, FermionField &src, FermionField &res,int dag);
#include <chroma.h>
#include <actions/ferm/invert/syssolver_linop_cg_array.h>
#include <actions/ferm/invert/syssolver_linop_aggregate.h>
namespace Chroma {
class ChromaWrapper {
public:
typedef multi1d<LatticeColorMatrix> U;
typedef LatticeStaggeredFermion T4;
static void ImportGauge(GaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd= gr._grid->GlobalDimensions();
for (x[0]=0;x[0]<gd[0];x[0]++){
for (x[1]=0;x[1]<gd[1];x[1]++){
for (x[2]=0;x[2]<gd[2];x[2]++){
for (x[3]=0;x[3]<gd[3];x[3]++){
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(LCM,gr,x);
for(int mu=0;mu<4;mu++){
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
cc = LCM(mu)()(i,j);
c = QDP::cmplx(QDP::Real(real(cc)),QDP::Real(imag(cc)));
QDP::pokeColor(cm,c,i,j);
}}
QDP::pokeSite(ch[mu],cm,cx);
}
}}}}
}
static void ExportGauge(GaugeField & gr,
QDP::multi1d<QDP::LatticeColorMatrix> & ch)
{
Grid::QCD::LorentzColourMatrix LCM;
Grid::Complex cc;
QDP::ColorMatrix cm;
QDP::Complex c;
std::vector<int> x(4);
QDP::multi1d<int> cx(4);
std::vector<int> gd= gr._grid->GlobalDimensions();
for (x[0]=0;x[0]<gd[0];x[0]++){
for (x[1]=0;x[1]<gd[1];x[1]++){
for (x[2]=0;x[2]<gd[2];x[2]++){
for (x[3]=0;x[3]<gd[3];x[3]++){
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
for(int mu=0;mu<4;mu++){
for(int i=0;i<3;i++){
for(int j=0;j<3;j++){
cm = QDP::peekSite(ch[mu],cx);
c = QDP::peekColor(cm,i,j);
cc = Grid::Complex(toDouble(real(c)),toDouble(imag(c)));
LCM(mu)()(i,j)= cc;
}}
}
Grid::pokeSite(LCM,gr,x);
}}}}
}
static void ImportFermion(FermionField & gr,
QDP::LatticeStaggeredFermion & ch )
{
Grid::QCD::ColourVector F;
Grid::Complex c;
std::vector<int> x(5);
QDP::multi1d<int> cx(4);
std::vector<int> gd= gr._grid->GlobalDimensions();
for (x[0]=0;x[0]<gd[0];x[0]++){
for (x[1]=0;x[1]<gd[1];x[1]++){
for (x[2]=0;x[2]<gd[2];x[2]++){
for (x[3]=0;x[3]<gd[3];x[3]++){
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
Grid::peekSite(F,gr,x);
QDP::ColorVector cv;
for(int j=0;j<3;j++){
QDP::Complex cc;
c = F()()(j) ;
cc = QDP::cmplx(QDP::Real(real(c)),QDP::Real(imag(c)));
pokeColor(cv,cc,j);
}
QDP::StaggeredFermion cF;
pokeSpin(cF,cv,0);
QDP::pokeSite(ch,cF,cx);
}}}}
}
static void ExportFermion(FermionField & gr,
QDP::LatticeStaggeredFermion & ch )
{
Grid::QCD::ColourVector F;
Grid::Complex c;
std::vector<int> x(5);
QDP::multi1d<int> cx(4);
std::vector<int> gd= gr._grid->GlobalDimensions();
for (x[0]=0;x[0]<gd[0];x[0]++){
for (x[1]=0;x[1]<gd[1];x[1]++){
for (x[2]=0;x[2]<gd[2];x[2]++){
for (x[3]=0;x[3]<gd[3];x[3]++){
cx[0] = x[0];
cx[1] = x[1];
cx[2] = x[2];
cx[3] = x[3];
QDP::StaggeredFermion cF = QDP::peekSite(ch,cx);
for(int j=0;j<3;j++){
QDP::ColorVector cS=QDP::peekSpin(cF,0);
QDP::Complex cc=QDP::peekColor(cS,j);
c = Grid::Complex(QDP::toDouble(QDP::real(cc)),
QDP::toDouble(QDP::imag(cc)));
F()()(j) = c;
}
Grid::pokeSite(F,gr,x);
}}}}
}
static Handle< Chroma::EvenOddLinearOperator<T4,U,U> > GetLinOp (U &u,U &u_fat,U &u_triple)
{
QDP::Real _mq(mq);
QDP::multi1d<int> bcs(QDP::Nd);
bcs[0] = bcs[1] = bcs[2] = bcs[3] = 1;
Chroma::AsqtadFermActParams p;
p.Mass = _mq;
p.u0 = Real(1.0);
Chroma::Handle<Chroma::FermBC<T4,U,U> > fbc(new Chroma::SimpleFermBC< T4, U, U >(bcs));
Chroma::Handle<Chroma::CreateFermState<T4,U,U> > cfs( new Chroma::CreateSimpleFermState<T4,U,U>(fbc));
Chroma::AsqtadFermAct S_f(cfs,p);
Chroma::Handle< Chroma::FermState<T4,U,U> > ffs( S_f.createState(u) );
u_fat =ffs.cast<AsqtadConnectStateBase>()->getFatLinks();
u_triple=ffs.cast<AsqtadConnectStateBase>()->getTripleLinks();
return S_f.linOp(ffs);
}
};
}
int main (int argc,char **argv )
{
/********************************************************
* Setup QDP
*********************************************************/
Chroma::initialize(&argc,&argv);
Chroma::WilsonTypeFermActs4DEnv::registerAll();
/********************************************************
* Setup Grid
*********************************************************/
Grid::Grid_init(&argc,&argv);
Grid::GridCartesian * UGrid = Grid::QCD::SpaceTimeGrid::makeFourDimGrid(Grid::GridDefaultLatt(),
Grid::GridDefaultSimd(Grid::QCD::Nd,Grid::vComplex::Nsimd()),
Grid::GridDefaultMpi());
std::vector<int> gd = UGrid->GlobalDimensions();
QDP::multi1d<int> nrow(QDP::Nd);
for(int mu=0;mu<4;mu++) nrow[mu] = gd[mu];
QDP::Layout::setLattSize(nrow);
QDP::Layout::create();
GaugeField uthin (UGrid);
GaugeField ufat (UGrid);
GaugeField utriple(UGrid);
FermionField src(UGrid);
FermionField res_chroma(UGrid);
FermionField res_grid (UGrid);
{
std::cout << "*****************************"<<std::endl;
std::cout << "Staggered Action " <<std::endl;
std::cout << "*****************************"<<std::endl;
make_gauge(uthin,src);
for(int dag=0;dag<2;dag++) {
std::cout << "Dag = "<<dag<<std::endl;
calc_chroma(uthin,utriple,ufat,src,res_chroma,dag);
// Remove the normalisation of Chroma Gauge links ??
std::cout << "Norm of chroma Asqtad multiply "<<Grid::norm2(res_chroma)<<std::endl;
calc_grid (uthin,utriple,ufat,src,res_grid,dag);
std::cout << "Norm of thin gauge "<< Grid::norm2(uthin) <<std::endl;
std::cout << "Norm of fat gauge "<< Grid::norm2(ufat) <<std::endl;
std::cout << "Norm of Grid Asqtad multiply "<<Grid::norm2(res_grid)<<std::endl;
/*
std::cout << " site 0 of Uthin "<<uthin._odata[0] <<std::endl;
std::cout << " site 0 of Utriple"<<utriple._odata[0] <<std::endl;
std::cout << " site 0 of Ufat "<<ufat._odata[0] <<std::endl;
std::cout << " site 0 of Grid "<<res_grid._odata[0] <<std::endl;
std::cout << " site 0 of Chroma "<<res_chroma._odata[0] <<std::endl;
*/
res_chroma=res_chroma - res_grid;
std::cout << "Norm of difference "<<Grid::norm2(res_chroma)<<std::endl;
}
}
std::cout << "Finished test "<<std::endl;
Chroma::finalize();
}
void calc_chroma(GaugeField & lat, GaugeField &uthin, GaugeField &ufat, FermionField &src, FermionField &res,int dag)
{
typedef QDP::LatticeStaggeredFermion T;
typedef QDP::multi1d<QDP::LatticeColorMatrix> U;
U u(4);
U ut(4);
U uf(4);
// Chroma::HotSt(u);
Chroma::ChromaWrapper::ImportGauge(lat,u) ;
QDP::LatticeStaggeredFermion check;
QDP::LatticeStaggeredFermion result;
QDP::LatticeStaggeredFermion tmp;
QDP::LatticeStaggeredFermion psi;
Chroma::ChromaWrapper::ImportFermion(src,psi);
auto linop =Chroma::ChromaWrapper::GetLinOp(u,uf,ut);
Chroma::ChromaWrapper::ExportGauge(uthin,ut) ;
Chroma::ChromaWrapper::ExportGauge(ufat ,uf) ;
enum Chroma::PlusMinus isign;
if ( dag ) {
isign=Chroma::MINUS;
} else {
isign=Chroma::PLUS;
}
std::cout << "Calling Chroma Linop "<< std::endl;
linop->evenEvenLinOp(tmp,psi,isign); check[rb[0]] = tmp;
linop->oddOddLinOp (tmp,psi,isign); check[rb[1]] = tmp;
linop->evenOddLinOp(tmp,psi,isign) ; check[rb[0]]+= tmp;
linop->oddEvenLinOp(tmp,psi,isign) ; check[rb[1]]+= tmp;
Chroma::ChromaWrapper::ExportFermion(res,check) ;
}
void make_gauge(GaugeField & Umu,FermionField &src)
{
using namespace Grid;
using namespace Grid::QCD;
std::vector<int> seeds4({1,2,3,4});
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Umu._grid;
Grid::GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
Grid::QCD::SU3::HotConfiguration(RNG4,Umu);
Grid::gaussian(RNG4,src);
}
void calc_grid(GaugeField & Uthin, GaugeField & Utriple, GaugeField & Ufat, FermionField &src, FermionField &res,int dag)
{
using namespace Grid;
using namespace Grid::QCD;
Grid::GridCartesian * UGrid = (Grid::GridCartesian *) Uthin._grid;
Grid::GridRedBlackCartesian * UrbGrid = Grid::QCD::SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
Grid::QCD::ImprovedStaggeredFermionR Dstag(Uthin,Utriple,Ufat,*UGrid,*UrbGrid,mq*2.0);
std::cout << Grid::GridLogMessage <<" Calling Grid staggered multiply "<<std::endl;
if ( dag )
Dstag.Mdag(src,res);
else
Dstag.M(src,res);
res = res ; // Convention mismatch to Chroma
return;
}

2
tests/solver/Test_cf_cr_unprec.cc
View File

@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_contfrac_cg.cc
View File

@ -94,7 +94,7 @@ int main (int argc, char ** argv)
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
RealD mass=0.1;

2
tests/solver/Test_dwf_cg_schur.cc
View File

@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_dwf_cg_unprec.cc
View File

@ -61,7 +61,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_dwf_cr_unprec.cc
View File

@ -65,7 +65,7 @@ int main (int argc, char ** argv)
LatticeFermion src(FGrid); random(RNG5,src);
LatticeFermion result(FGrid); result=zero;
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(RNG4,Umu);
std::vector<LatticeColourMatrix> U(4,UGrid);

2
tests/solver/Test_wilson_cg_prec.cc
View File

@ -60,7 +60,7 @@ int main (int argc, char ** argv)
LatticeFermion src(&Grid); random(pRNG,src);
RealD nrm = norm2(src);
LatticeFermion result(&Grid); result=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);

2
tests/solver/Test_wilson_cg_schur.cc
View File

@ -57,7 +57,7 @@ int main (int argc, char ** argv)
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds);
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
LatticeFermion src(&Grid); random(pRNG,src);
LatticeFermion result(&Grid); result=zero;

2
tests/solver/Test_wilson_cg_unprec.cc
View File

@ -60,7 +60,7 @@ int main (int argc, char ** argv)
LatticeFermion src(&Grid); random(pRNG,src);
RealD nrm = norm2(src);
LatticeFermion result(&Grid); result=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
double volume=1;
for(int mu=0;mu<Nd;mu++){

2
tests/solver/Test_wilson_cr_unprec.cc
View File

@ -60,7 +60,7 @@ int main (int argc, char ** argv)
LatticeFermion src(&Grid); random(pRNG,src);
RealD nrm = norm2(src);
LatticeFermion result(&Grid); result=zero;
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
std::vector<LatticeColourMatrix> U(4,&Grid);

113
tests/solver/Test_zmobius_cg_prec.cc Normal file
View File

@ -0,0 +1,113 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_cg_prec.cc
Copyright (C) 2015
Author: Peter Boyle <paboyle@ph.ed.ac.uk>
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License along
with this program; if not, write to the Free Software Foundation, Inc.,
51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
See the full license in the file "LICENSE" in the top level distribution
directory
*************************************************************************************/
/* END LEGAL */
#include <Grid/Grid.h>
using namespace std;
using namespace Grid;
using namespace Grid::QCD;
template <class d>
struct scal {
d internal;
};
Gamma::Algebra Gmu[] = {Gamma::Algebra::GammaX, Gamma::Algebra::GammaY, Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT};
int main(int argc, char** argv) {
Grid_init(&argc, &argv);
const int Ls = 16;
GridCartesian* UGrid = SpaceTimeGrid::makeFourDimGrid(
GridDefaultLatt(), GridDefaultSimd(Nd, vComplex::Nsimd()),
GridDefaultMpi());
GridRedBlackCartesian* UrbGrid =
SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
GridCartesian* FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls, UGrid);
GridRedBlackCartesian* FrbGrid =
SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls, UGrid);
std::vector<int> seeds4({1, 2, 3, 4});
std::vector<int> seeds5({5, 6, 7, 8});
GridParallelRNG RNG5(FGrid);
RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid);
RNG4.SeedFixedIntegers(seeds4);
LatticeFermion src(FGrid);
random(RNG5, src);
LatticeFermion result(FGrid);
result = zero;
LatticeGaugeField Umu(UGrid);
SU3::HotConfiguration(RNG4, Umu);
std::cout << GridLogMessage << "Lattice dimensions: " << GridDefaultLatt()
<< " Ls: " << Ls << std::endl;
std::vector<LatticeColourMatrix> U(4, UGrid);
for (int mu = 0; mu < Nd; mu++) {
U[mu] = PeekIndex<LorentzIndex>(Umu, mu);
}
RealD mass = 0.01;
RealD M5 = 1.8;
std::vector < std::complex<double> > omegas;
for(int i=0;i<Ls;i++){
double imag = 0.;
if (i==0) imag=1.;
if (i==Ls-1) imag=-1.;
std::complex<double> temp (0.25+0.01*i, imag*0.01);
omegas.push_back(temp);
}
ZMobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,1.,0.);
LatticeFermion src_o(FrbGrid);
LatticeFermion result_o(FrbGrid);
pickCheckerboard(Odd, src_o, src);
result_o = zero;
GridStopWatch CGTimer;
SchurDiagMooeeOperator<ZMobiusFermionR, LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8, 10000, 0);// switch off the assert
CGTimer.Start();
CG(HermOpEO, src_o, result_o);
CGTimer.Stop();
std::cout << GridLogMessage << "Total CG time : " << CGTimer.Elapsed()
<< std::endl;
std::cout << GridLogMessage << "######## Dhop calls summary" << std::endl;
Ddwf.Report();
Grid_finalize();
}