portelli/Grid
1
0
Fork 0
mirror of https://github.com/paboyle/Grid.git synced 2025-06-13 04:37:05 +01:00
Code Releases Activity

Merge branch 'develop' into feature/scalar_adjointFT

Browse Source
This commit is contained in:
Guido Cossu
2017-05-05 15:47:33 +01:00
parent 7b03d8d087 47b5c07ffb
commit 8546d01a4c
115 changed files with 7439 additions and 3818 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
tests/IO/Test_serialisation.cc
View File

@ -115,6 +115,7 @@ int main(int argc,char **argv)
// test serializable class writing
myclass obj(1234); // non-trivial constructor
std::vector<myclass> vec;
std::pair<myenum, myenum> pair;
std::cout << "-- serialisable class writing to 'bother.xml'..." << std::endl;
write(WR,"obj",obj);
@ -122,6 +123,8 @@ int main(int argc,char **argv)
vec.push_back(myclass(1234));
vec.push_back(myclass(5678));
vec.push_back(myclass(3838));
pair = std::make_pair(myenum::red, myenum::blue);
write(WR, "objvec", vec);
std::cout << "-- serialisable class writing to std::cout:" << std::endl;
std::cout << obj << std::endl;
@ -129,21 +132,30 @@ int main(int argc,char **argv)
std::cout << "vec[0] == obj: " << ((vec[0] == obj) ? "true" : "false") << std::endl;
std::cout << "vec[1] == obj: " << ((vec[1] == obj) ? "true" : "false") << std::endl;
write(WR, "objpair", pair);
std::cout << "-- pair writing to std::cout:" << std::endl;
std::cout << pair << std::endl;
// read tests
std::cout << "\n==== IO self-consistency tests" << std::endl;
//// XML
ioTest<XmlWriter, XmlReader>("iotest.xml", obj, "XML (object) ");
ioTest<XmlWriter, XmlReader>("iotest.xml", vec, "XML (vector of objects)");
ioTest<XmlWriter, XmlReader>("iotest.xml", pair, "XML (pair of objects)");
//// binary
ioTest<BinaryWriter, BinaryReader>("iotest.bin", obj, "binary (object) ");
ioTest<BinaryWriter, BinaryReader>("iotest.bin", vec, "binary (vector of objects)");
ioTest<BinaryWriter, BinaryReader>("iotest.bin", pair, "binary (pair of objects)");
//// text
ioTest<TextWriter, TextReader>("iotest.dat", obj, "text (object) ");
ioTest<TextWriter, TextReader>("iotest.dat", vec, "text (vector of objects)");
ioTest<TextWriter, TextReader>("iotest.dat", pair, "text (pair of objects)");
//// HDF5
#undef HAVE_HDF5
#ifdef HAVE_HDF5
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", obj, "HDF5 (object) ");
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", vec, "HDF5 (vector of objects)");
ioTest<Hdf5Writer, Hdf5Reader>("iotest.h5", pair, "HDF5 (pair of objects)");
#endif
std::cout << "\n==== vector flattening/reconstruction" << std::endl;

110
tests/Test_simd.cc
View File

@ -308,18 +308,23 @@ 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 {
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 ];
for(int i=0;i<sz;i++) {
int j1 = i&(~msk);
int j2 = i|msk;
if ( (i&msk) == 0 ) { r1[i]=in1[j1];}
else { r1[i]=in2[j1];}
if ( (i&msk) == 0 ) { r2[i]=in1[j2];}
else { r2[i]=in2[j2];}
}
}
std::string name(void) const { return std::string("Exchange"); }
};
@ -454,8 +459,8 @@ void ExchangeTester(const functor &func)
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++) std::cout << " i "<<i<<" ref "<<reference1[i]<<" res "<<result1[i]<<std::endl;
//for(int i=0;i<Nsimd;i++) std::cout << " i "<<i<<" ref "<<reference2[i]<<" res "<<result2[i]<<std::endl;
for(int i=0;i<Nsimd;i++){
int found=0;
@ -465,7 +470,7 @@ void ExchangeTester(const functor &func)
// std::cout << " i "<<i<<" j "<<j<<" "<<reference1[j]<<" "<<result1[i]<<std::endl;
}
}
assert(found==1);
// assert(found==1);
}
for(int i=0;i<Nsimd;i++){
int found=0;
@ -475,15 +480,24 @@ void ExchangeTester(const functor &func)
// std::cout << " i "<<i<<" j "<<j<<" "<<reference2[j]<<" "<<result2[i]<<std::endl;
}
}
assert(found==1);
// assert(found==1);
}
/*
for(int i=0;i<Nsimd;i++){
std::cout << " i "<< i
<<" result1 "<<result1[i]
<<" result2 "<<result2[i]
<<" test1 "<<test1[i]
<<" test2 "<<test2[i]
<<" input1 "<<input1[i]
<<" input2 "<<input2[i]<<std::endl;
}
*/
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;
// }
}
}
@ -678,5 +692,69 @@ int main (int argc, char ** argv)
IntTester(funcMinus());
IntTester(funcTimes());
std::cout<<GridLogMessage << "==================================="<< std::endl;
std::cout<<GridLogMessage << "Testing precisionChange "<< std::endl;
std::cout<<GridLogMessage << "==================================="<< std::endl;
{
GridSerialRNG sRNG;
sRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
const int Ndp = 16;
const int Nsp = Ndp/2;
const int Nhp = Ndp/4;
std::vector<vRealH,alignedAllocator<vRealH> > H (Nhp);
std::vector<vRealF,alignedAllocator<vRealF> > F (Nsp);
std::vector<vRealF,alignedAllocator<vRealF> > FF(Nsp);
std::vector<vRealD,alignedAllocator<vRealD> > D (Ndp);
std::vector<vRealD,alignedAllocator<vRealD> > DD(Ndp);
for(int i=0;i<16;i++){
random(sRNG,D[i]);
}
// Double to Single
precisionChange(&F[0],&D[0],Ndp);
precisionChange(&DD[0],&F[0],Ndp);
std::cout << GridLogMessage<<"Double to single";
for(int i=0;i<Ndp;i++){
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i] <<std::endl;
DD[i] = DD[i] - D[i];
decltype(innerProduct(DD[0],DD[0])) nrm;
nrm = innerProduct(DD[i],DD[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-14 );
}
std::cout <<" OK ! "<<std::endl;
// Double to Half
#ifdef USE_FP16
std::cout << GridLogMessage<< "Double to half" ;
precisionChange(&H[0],&D[0],Ndp);
precisionChange(&DD[0],&H[0],Ndp);
for(int i=0;i<Ndp;i++){
// std::cout << "DD["<<i<<"] = "<< DD[i]<<" "<<D[i]<<" "<<DD[i]-D[i]<<std::endl;
DD[i] = DD[i] - D[i];
decltype(innerProduct(DD[0],DD[0])) nrm;
nrm = innerProduct(DD[i],DD[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-3 );
}
std::cout <<" OK ! "<<std::endl;
std::cout << GridLogMessage<< "Single to half";
// Single to Half
precisionChange(&H[0] ,&F[0],Nsp);
precisionChange(&FF[0],&H[0],Nsp);
for(int i=0;i<Nsp;i++){
// std::cout << "FF["<<i<<"] = "<< FF[i]<<" "<<F[i]<<" "<<FF[i]-F[i]<<std::endl;
FF[i] = FF[i] - F[i];
decltype(innerProduct(FF[0],FF[0])) nrm;
nrm = innerProduct(FF[i],FF[i]);
auto tmp = Reduce(nrm);
// std::cout << tmp << std::endl;
assert( tmp < 1.0e-3 );
}
std::cout <<" OK ! "<<std::endl;
#endif
}
Grid_finalize();
}

4
tests/core/Test_fft_gfix.cc
View File

@ -148,11 +148,13 @@ class FourierAcceleratedGaugeFixer : public Gimpl {
Complex psqMax(16.0);
Fp = psqMax*one/psq;
/*
static int once;
if ( once == 0 ) {
std::cout << " Fp " << Fp <<std::endl;
once ++;
}
}*/
pokeSite(TComplex(1.0),Fp,coor);
dmuAmu_p = dmuAmu_p * Fp;

287
tests/core/Test_mobius_even_odd.cc Normal file
View File

@ -0,0 +1,287 @@
/*************************************************************************************
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=10;
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;
#if 0
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);
}
#else
omegas.push_back( std::complex<double>(1.45806438985048,-0) );
omegas.push_back( std::complex<double>(1.18231318389348,-0) );
omegas.push_back( std::complex<double>(0.830951166685955,-0) );
omegas.push_back( std::complex<double>(0.542352409156791,-0) );
omegas.push_back( std::complex<double>(0.341985020453729,-0) );
omegas.push_back( std::complex<double>(0.21137902619029,-0) );
omegas.push_back( std::complex<double>(0.126074299502912,-0) );
omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
#endif
MobiusFermionR Ddwf(Umu, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, 0.5,0.5);
// 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<MobiusFermionR,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();
}

9
tests/core/Test_wilson_even_odd.cc
View File

@ -2,11 +2,10 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_even_odd.cc
Source file: ./tests/Test_wilson_tm_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
@ -89,8 +88,8 @@ int main (int argc, char ** argv)
}
RealD mass=0.1;
RealD mu = 0.1;
WilsonTMFermionR Dw(Umu,Grid,RBGrid,mass,mu);
WilsonFermionR Dw(Umu,Grid,RBGrid,mass);
LatticeFermion src_e (&RBGrid);
LatticeFermion src_o (&RBGrid);
@ -207,7 +206,7 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
SchurDiagMooeeOperator<WilsonTMFermionR,LatticeFermion> HermOpEO(Dw);
SchurDiagMooeeOperator<WilsonFermionR,LatticeFermion> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);

9
tests/core/Test_wilson_tm_even_odd.cc → tests/core/Test_wilson_twisted_mass_even_odd.cc
View File

@ -2,10 +2,11 @@
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_tm_even_odd.cc
Source file: ./tests/Test_wilson_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
@ -88,8 +89,8 @@ int main (int argc, char ** argv)
}
RealD mass=0.1;
WilsonFermionR Dw(Umu,Grid,RBGrid,mass);
RealD mu = 0.1;
WilsonTMFermionR Dw(Umu,Grid,RBGrid,mass,mu);
LatticeFermion src_e (&RBGrid);
LatticeFermion src_o (&RBGrid);
@ -206,7 +207,7 @@ int main (int argc, char ** argv)
pickCheckerboard(Odd ,phi_o,phi);
RealD t1,t2;
SchurDiagMooeeOperator<WilsonFermionR,LatticeFermion> HermOpEO(Dw);
SchurDiagMooeeOperator<WilsonTMFermionR,LatticeFermion> HermOpEO(Dw);
HermOpEO.MpcDagMpc(chi_e,dchi_e,t1,t2);
HermOpEO.MpcDagMpc(chi_o,dchi_o,t1,t2);

13
tests/debug/Test_synthetic_lanczos.cc
View File

@ -115,8 +115,8 @@ int main (int argc, char ** argv)
RNG.SeedFixedIntegers(seeds);
RealD alpha = 1.0;
RealD beta = 0.03;
RealD alpha = 1.2;
RealD beta = 0.1;
RealD mu = 0.0;
int order = 11;
ChebyshevLanczos<LatticeComplex> Cheby(alpha,beta,mu,order);
@ -131,10 +131,9 @@ int main (int argc, char ** argv)
const int Nit= 10000;
int Nconv;
RealD eresid = 1.0e-8;
RealD eresid = 1.0e-6;
ImplicitlyRestartedLanczos<LatticeComplex> IRL(HermOp,X,Nk,Nm,eresid,Nit);
ImplicitlyRestartedLanczos<LatticeComplex> ChebyIRL(HermOp,Cheby,Nk,Nm,eresid,Nit);
LatticeComplex src(grid); gaussian(RNG,src);
@ -145,9 +144,9 @@ int main (int argc, char ** argv)
}
{
// std::vector<RealD> eval(Nm);
// std::vector<LatticeComplex> evec(Nm,grid);
// ChebyIRL.calc(eval,evec,src, Nconv);
std::vector<RealD> eval(Nm);
std::vector<LatticeComplex> evec(Nm,grid);
ChebyIRL.calc(eval,evec,src, Nconv);
}
Grid_finalize();

15
tests/forces/Test_dwf_force.cc
View File

@ -45,8 +45,19 @@ int main (int argc, char ** argv)
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});
// Want a different conf at every run
// First create an instance of an engine.
std::random_device rnd_device;
// Specify the engine and distribution.
std::mt19937 mersenne_engine(rnd_device());
std::uniform_int_distribution<int> dist(1, 100);
auto gen = std::bind(dist, mersenne_engine);
std::vector<int> seeds4(4);
generate(begin(seeds4), end(seeds4), gen);
//std::vector<int> seeds4({1,2,3,5});
std::vector<int> seeds5({5,6,7,8});
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);

153
tests/forces/Test_mobius_force.cc Normal file
View File

@ -0,0 +1,153 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_force.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;
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=8;
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);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
LatticeFermion phi (FGrid); gaussian(RNG5,phi);
LatticeFermion Mphi (FGrid);
LatticeFermion MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD mass=0.01;
RealD M5=1.8;
RealD b=0.5;
RealD c=0.5;
MobiusFermionR Ddwf(U,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,b,c);
Ddwf.M (phi,Mphi);
ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(UGrid);
LatticeGaugeField tmp(UGrid);
Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
LatticeFermion Ftmp (FGrid);
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
// fourth order exponential approx
parallel_for(auto i=mom.begin();i<mom.end();i++){
Uprime[i](mu) =
U[i](mu)
+ mom[i](mu)*U[i](mu)*dt
+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
;
}
}
Ddwf.ImportGauge(Uprime);
Ddwf.M (phi,MphiPrime);
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(UGrid); dS = zero;
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
}
Complex dSpred = sum(dS);
std::cout << GridLogMessage << " -- S "<<S<<std::endl;
std::cout << GridLogMessage << " -- Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

167
tests/forces/Test_wilson_force_phiMdagMphi.cc
View File

@ -1,167 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_force_phiMdagMphi.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;
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::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
LatticeFermion phi (&Grid); gaussian(pRNG,phi);
LatticeFermion Mphi (&Grid);
LatticeFermion Mdagphi (&Grid);
LatticeFermion MphiPrime (&Grid);
LatticeFermion MdagphiPrime (&Grid);
LatticeFermion dMphi (&Grid);
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
// SU3::ColdConfiguration(pRNG,U);
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD mass=-4.0; //kills the diagonal term
WilsonFermionR Dw (U, Grid,RBGrid,mass);
Dw.M (phi,Mphi);
Dw.Mdag(phi,Mdagphi);
ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p
ComplexD Sdag = innerProduct(Mdagphi,Mdagphi); // pdag MMdag p
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(&Grid);
LatticeGaugeField UdSdUdag(&Grid);
LatticeGaugeField tmp(&Grid);
Dw.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
Dw.MDeriv(tmp , Mdagphi, phi,DaggerYes ); UdSdUdag=tmp;
LatticeFermion dMdagphi (&Grid); dMdagphi=zero;
LatticeFermion Ftmp (&Grid);
// Dw.MDeriv(UdSdU,Mdagphi, phi,DaggerYes );// UdSdU =UdSdU +tmp;
////////////////////////////////////
// Modify the gauge field a little in one dir
////////////////////////////////////
RealD dt = 1.0e-3;
LatticeColourMatrix mommu(&Grid);
LatticeGaugeField mom(&Grid);
LatticeGaugeField Uprime(&Grid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mommu); // Traceless antihermitian momentum; gaussian in lie alg
// Dw.DoubleStore(Dw.Umu,Uprime); // update U _and_ Udag
Dw.DhopDirDisp(phi,Ftmp,mu,mu+4,DaggerYes);
dMdagphi=dMdagphi+mommu*Ftmp*dt;
PokeIndex<LorentzIndex>(mom,mommu,mu);
parallel_for(auto i=mom.begin();i<mom.end();i++){
Uprime[i](mu) =U[i](mu)+ mom[i](mu)*U[i](mu)*dt;
Dw.Umu[i](mu) =Uprime[i](mu); // update U but _not_ Udag
}
}
Dw.Mdag(phi,MdagphiPrime);
Dw.M (phi,MphiPrime);
std::cout << GridLogMessage << "deltaMdag phi "<< norm2(dMdagphi) <<std::endl;
Ftmp=MdagphiPrime - Mdagphi;
std::cout << GridLogMessage << "diff Mdag phi "<< norm2(Ftmp) <<std::endl;
Ftmp = Ftmp - dMdagphi;
std::cout << GridLogMessage << "err Mdag phi "<< norm2(Ftmp) <<std::endl;
std::cout << dMdagphi<<std::endl;
Ftmp=MdagphiPrime - Mdagphi;
std::cout << Ftmp<<std::endl;
ComplexD Sprime = innerProduct(Mphi ,MphiPrime);
ComplexD Sprimedag = innerProduct(Mdagphi,MdagphiPrime);
ComplexD deltaSdag = innerProduct(Mdagphi,dMdagphi);
std::cout << GridLogMessage << "deltaSdag from inner prod of mom* M[u] "<<deltaSdag<<std::endl;
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(&Grid); dS = zero;
LatticeComplex dSdag(&Grid); dSdag = zero;
parallel_for(auto i=mom.begin();i<mom.end();i++){
for(int mu=0;mu<Nd;mu++){
// dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) - mom[i](mu)* adj( UdSdU[i](mu)) )*dt;
dS[i]() = dS[i]()+trace(mom[i](mu) * UdSdU[i](mu) )*dt;
dSdag[i]() = dSdag[i]()+trace(mom[i](mu) * UdSdUdag[i](mu) )*dt;
}
}
Complex dSpred = sum(dS);
Complex dSdagpred = sum(dSdag);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout << "\n\n"<<std::endl;
std::cout << GridLogMessage << " Sdag "<<Sdag<<std::endl;
std::cout << GridLogMessage << " Sprimedag "<<Sprimedag<<std::endl;
std::cout << GridLogMessage << "dSdag "<<Sprimedag-Sdag<<std::endl;
std::cout << GridLogMessage << "predict dSdag "<< dSdagpred <<std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

189
tests/forces/Test_wilson_force_phiMphi.cc
View File

@ -1,189 +0,0 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_force_phiMphi.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;
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::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid);
pRNG.SeedFixedIntegers(std::vector<int>({45,12,81,9}));
LatticeFermion phi (&Grid); gaussian(pRNG,phi);
LatticeFermion Mphi (&Grid);
LatticeFermion MphiPrime (&Grid);
LatticeFermion dMphi (&Grid);
LatticeGaugeField U(&Grid);
SU3::HotConfiguration(pRNG,U);
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD mass=-4.0; //kills the diagonal term
WilsonFermionR Dw (U, Grid,RBGrid,mass);
Dw.M(phi,Mphi);
ComplexD S = innerProduct(phi,Mphi);
// get the deriv
LatticeGaugeField UdSdU(&Grid);
Dw.MDeriv(UdSdU,phi, phi,DaggerNo );
////////////////////////////////////
// Modify the gauge field a little in one dir
////////////////////////////////////
RealD dt = 1.0e-3;
Complex Complex_i(0,1);
LatticeColourMatrix Umu(&Grid);
LatticeColourMatrix Umu_save(&Grid);
LatticeColourMatrix dU (&Grid);
LatticeColourMatrix mom(&Grid);
SU3::GaussianFundamentalLieAlgebraMatrix(pRNG, mom); // Traceless antihermitian momentum; gaussian in lie alg
// check mom is as i expect
LatticeColourMatrix tmpmom(&Grid);
tmpmom = mom+adj(mom);
std::cout << GridLogMessage << "mom anti-herm check "<< norm2(tmpmom)<<std::endl;
std::cout << GridLogMessage << "mom tr check "<< norm2(trace(mom))<<std::endl;
const int mu=0;
Umu = PeekIndex<LorentzIndex>(U,mu);
Umu_save=Umu;
dU = mom * Umu * dt;
Umu= Umu+dU;
PokeIndex<LorentzIndex>(Dw.Umu,Umu,mu);
Dw.M(phi,MphiPrime);
ComplexD Sprime = innerProduct(phi,MphiPrime);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
Dw.Umu=zero;
PokeIndex<LorentzIndex>(Dw.Umu,dU,mu);
Dw.M(phi,dMphi);
ComplexD deltaS = innerProduct(phi,dMphi);
std::cout << GridLogMessage << "deltaS "<<deltaS<<std::endl;
Dw.Umu=zero;
PokeIndex<LorentzIndex>(Dw.Umu,Umu_save,mu);
Dw.Mdir(phi,dMphi,mu,1);
dMphi = dt*mom*dMphi;
deltaS = innerProduct(phi,dMphi);
std::cout << GridLogMessage << "deltaS from inner prod of mom* M[u] "<<deltaS<<std::endl;
deltaS = sum(trace(outerProduct(dMphi,phi)));
std::cout << GridLogMessage << "deltaS from trace outer prod of deltaM "<<deltaS<<std::endl;
/*
LatticeComplex lip(&Grid);
lip = localInnerProduct(phi,dMphi);
LatticeComplex trop(&Grid);
trop = trace(outerProduct(dMphi,phi));
LatticeSpinColourMatrix op(&Grid);
op = outerProduct(dMphi,phi);
LatticeSpinColourMatrix hop(&Grid);
LatticeComplex op_cpt(&Grid);
for(int s1=0;s1<Ns;s1++){
for(int s2=0;s2<Ns;s2++){
for(int c1=0;c1<Nc;c1++){
for(int c2=0;c2<Nc;c2++){
op_cpt = peekColour(peekSpin(dMphi,s1),c1) * adj(peekColour(peekSpin(phi,s2),c2));
parallel_for(auto i=hop.begin();i<hop.end();i++){
hop[i]()(s1,s2)(c1,c2) = op_cpt[i]()()();
}
}}}}
LatticeSpinColourMatrix diffop(&Grid);
diffop = hop - op;
std::cout << GridLogMessage << "hand outer prod diff "<<norm2(diffop)<<std::endl;
deltaS = sum(trace(hop));
std::cout << GridLogMessage << "deltaS hop "<<deltaS<<std::endl;
std::cout << GridLogMessage<< " phi[0] : "<< phi._odata[0]<<std::endl;
std::cout << GridLogMessage<< "dMphi[0] : "<<dMphi._odata[0]<<std::endl;
std::cout << GridLogMessage<< "hop[0] : "<< hop._odata[0]<<std::endl;
std::cout << GridLogMessage<< " op[0] : "<< op._odata[0]<<std::endl;
std::cout << GridLogMessage << "lip "<<lip<<std::endl;
std::cout << GridLogMessage << "trop "<<trop<<std::endl;
*/
// std::cout << GridLogMessage << " UdSdU " << UdSdU << std::endl;
LatticeComplex dS(&Grid); dS = zero;
parallel_for(auto i=mom.begin();i<mom.end();i++){
dS[i]() = trace(mom[i]() * UdSdU[i](mu) )*dt;
}
Complex dSpred = sum(dS);
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

167
tests/forces/Test_zmobius_force.cc Normal file
View File

@ -0,0 +1,167 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_dwf_force.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;
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=8;
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);
int threads = GridThread::GetThreads();
std::cout<<GridLogMessage << "Grid is setup to use "<<threads<<" threads"<<std::endl;
LatticeFermion phi (FGrid); gaussian(RNG5,phi);
LatticeFermion Mphi (FGrid);
LatticeFermion MphiPrime (FGrid);
LatticeGaugeField U(UGrid);
SU3::HotConfiguration(RNG4,U);
////////////////////////////////////
// Unmodified matrix element
////////////////////////////////////
RealD mass=0.01;
RealD M5=1.8;
RealD b=0.5;
RealD c=0.5;
std::vector < std::complex<double> > omegas;
omegas.push_back( std::complex<double>(1.45806438985048,-0) );
omegas.push_back( std::complex<double>(1.18231318389348,-0) );
omegas.push_back( std::complex<double>(0.830951166685955,-0) );
omegas.push_back( std::complex<double>(0.542352409156791,-0) );
omegas.push_back( std::complex<double>(0.341985020453729,-0) );
omegas.push_back( std::complex<double>(0.21137902619029,-0) );
omegas.push_back( std::complex<double>(0.126074299502912,-0) );
omegas.push_back( std::complex<double>(0.0990136651962626,-0) );
omegas.push_back( std::complex<double>(0.0686324988446592,0.0550658530827402) );
omegas.push_back( std::complex<double>(0.0686324988446592,-0.0550658530827402) );
ZMobiusFermionR Ddwf(U, *FGrid, *FrbGrid, *UGrid, *UrbGrid, mass, M5, omegas,b,c);
Ddwf.M (phi,Mphi);
ComplexD S = innerProduct(Mphi,Mphi); // pdag MdagM p
// get the deriv of phidag MdagM phi with respect to "U"
LatticeGaugeField UdSdU(UGrid);
LatticeGaugeField tmp(UGrid);
Ddwf.MDeriv(tmp , Mphi, phi,DaggerNo ); UdSdU=tmp;
Ddwf.MDeriv(tmp , phi, Mphi,DaggerYes ); UdSdU=(UdSdU+tmp);
LatticeFermion Ftmp (FGrid);
////////////////////////////////////
// Modify the gauge field a little
////////////////////////////////////
RealD dt = 0.0001;
LatticeColourMatrix mommu(UGrid);
LatticeColourMatrix forcemu(UGrid);
LatticeGaugeField mom(UGrid);
LatticeGaugeField Uprime(UGrid);
for(int mu=0;mu<Nd;mu++){
SU3::GaussianFundamentalLieAlgebraMatrix(RNG4, mommu); // Traceless antihermitian momentum; gaussian in lie alg
PokeIndex<LorentzIndex>(mom,mommu,mu);
// fourth order exponential approx
parallel_for(auto i=mom.begin();i<mom.end();i++){
Uprime[i](mu) =
U[i](mu)
+ mom[i](mu)*U[i](mu)*dt
+ mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt/2.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt/6.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt/24.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt/120.0)
+ mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *mom[i](mu) *U[i](mu)*(dt*dt*dt*dt*dt*dt/720.0)
;
}
}
Ddwf.ImportGauge(Uprime);
Ddwf.M (phi,MphiPrime);
ComplexD Sprime = innerProduct(MphiPrime ,MphiPrime);
//////////////////////////////////////////////
// Use derivative to estimate dS
//////////////////////////////////////////////
LatticeComplex dS(UGrid); dS = zero;
for(int mu=0;mu<Nd;mu++){
mommu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu=Ta(mommu)*2.0;
PokeIndex<LorentzIndex>(UdSdU,mommu,mu);
}
for(int mu=0;mu<Nd;mu++){
forcemu = PeekIndex<LorentzIndex>(UdSdU,mu);
mommu = PeekIndex<LorentzIndex>(mom,mu);
// Update PF action density
dS = dS+trace(mommu*forcemu)*dt;
}
Complex dSpred = sum(dS);
std::cout << GridLogMessage << " S "<<S<<std::endl;
std::cout << GridLogMessage << " Sprime "<<Sprime<<std::endl;
std::cout << GridLogMessage << "dS "<<Sprime-S<<std::endl;
std::cout << GridLogMessage << "predict dS "<< dSpred <<std::endl;
std::cout<< GridLogMessage << "Done" <<std::endl;
Grid_finalize();
}

368
tests/hadrons/Test_hadrons.hpp Normal file
View File

@ -0,0 +1,368 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons.hpp
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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.
*******************************************************************************/
#include <Grid/Hadrons/Application.hpp>
using namespace Grid;
using namespace Hadrons;
/*******************************************************************************
* Macros to reduce code duplication.
******************************************************************************/
// Useful definitions
#define ZERO_MOM "0. 0. 0. 0."
#define INIT_INDEX(s, n) (std::string(s) + "_" + std::to_string(n))
#define ADD_INDEX(s, n) (s + "_" + std::to_string(n))
#define LABEL_3PT(s, t1, t2) ADD_INDEX(INIT_INDEX(s, t1), t2)
#define LABEL_4PT(s, t1, t2, t3) ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3)
#define LABEL_4PT_NOISE(s, t1, t2, t3, nn) ADD_INDEX(ADD_INDEX(ADD_INDEX(INIT_INDEX(s, t1), t2), t3), nn)
// Wall source/sink macros
#define NAME_3MOM_WALL_SOURCE(t, mom) ("wall_" + std::to_string(t) + "_" + mom)
#define NAME_WALL_SOURCE(t) NAME_3MOM_WALL_SOURCE(t, ZERO_MOM)
#define NAME_POINT_SOURCE(pos) ("point_" + pos)
#define MAKE_3MOM_WALL_PROP(tW, mom, propName, solver)\
{\
std::string srcName = NAME_3MOM_WALL_SOURCE(tW, mom);\
makeWallSource(application, srcName, tW, mom);\
makePropagator(application, propName, srcName, solver);\
}
#define MAKE_WALL_PROP(tW, propName, solver)\
MAKE_3MOM_WALL_PROP(tW, ZERO_MOM, propName, solver)
// Sequential source macros
#define MAKE_SEQUENTIAL_PROP(tS, qSrc, mom, propName, solver)\
{\
std::string srcName = ADD_INDEX(qSrc + "_seq", tS);\
makeSequentialSource(application, srcName, qSrc, tS, mom);\
makePropagator(application, propName, srcName, solver);\
}
// Point source macros
#define MAKE_POINT_PROP(pos, propName, solver)\
{\
std::string srcName = NAME_POINT_SOURCE(pos);\
makePointSource(application, srcName, pos);\
makePropagator(application, propName, srcName, solver);\
}
/*******************************************************************************
* Functions for propagator construction.
******************************************************************************/
/*******************************************************************************
* Name: makePointSource
* Purpose: Construct point source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* pos - Position of point source.
* Returns: None.
******************************************************************************/
inline void makePointSource(Application &application, std::string srcName,
std::string pos)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::Point::Par pointPar;
pointPar.position = pos;
application.createModule<MSource::Point>(srcName, pointPar);
}
}
/*******************************************************************************
* Name: makeSequentialSource
* Purpose: Construct sequential source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* qSrc - Input quark for sequential inversion.
* tS - sequential source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeSequentialSource(Application &application, std::string srcName,
std::string qSrc, unsigned int tS,
std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::SeqGamma::Par seqPar;
seqPar.q = qSrc;
seqPar.tA = tS;
seqPar.tB = tS;
seqPar.mom = mom;
seqPar.gamma = Gamma::Algebra::GammaT;
application.createModule<MSource::SeqGamma>(srcName, seqPar);
}
}
/*******************************************************************************
* Name: makeWallSource
* Purpose: Construct wall source and add to application module.
* Parameters: application - main application that stores modules.
* srcName - name of source module to create.
* tW - wall source timeslice.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeWallSource(Application &application, std::string srcName,
unsigned int tW, std::string mom = ZERO_MOM)
{
// If the source already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(srcName)))
{
MSource::Wall::Par wallPar;
wallPar.tW = tW;
wallPar.mom = mom;
application.createModule<MSource::Wall>(srcName, wallPar);
}
}
/*******************************************************************************
* Name: makeWallSink
* Purpose: Wall sink smearing of a propagator.
* Parameters: application - main application that stores modules.
* propName - name of input propagator.
* wallName - name of smeared propagator.
* mom - momentum insertion (default is zero).
* Returns: None.
******************************************************************************/
inline void makeWallSink(Application &application, std::string propName,
std::string wallName, std::string mom = ZERO_MOM)
{
// If the propagator has already been smeared, don't smear it again.
// Temporarily removed, strategy for sink smearing likely to change.
/*if (!(Environment::getInstance().hasModule(wallName)))
{
MSink::Wall::Par wallPar;
wallPar.q = propName;
wallPar.mom = mom;
application.createModule<MSink::Wall>(wallName, wallPar);
}*/
}
/*******************************************************************************
* Name: makePropagator
* Purpose: Construct source and propagator then add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of source module to use.
* solver - solver to use (default is CG).
* Returns: None.
******************************************************************************/
inline void makePropagator(Application &application, std::string &propName,
std::string &srcName, std::string &solver)
{
// If the propagator already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(propName)))
{
Quark::Par quarkPar;
quarkPar.source = srcName;
quarkPar.solver = solver;
application.createModule<Quark>(propName, quarkPar);
}
}
/*******************************************************************************
* Name: makeLoop
* Purpose: Use noise source and inversion result to make loop propagator, then
* add to application module.
* Parameters: application - main application that stores modules.
* propName - name of propagator module to create.
* srcName - name of noise source module to use.
* resName - name of inversion result on given noise source.
* Returns: None.
******************************************************************************/
inline void makeLoop(Application &application, std::string &propName,
std::string &srcName, std::string &resName)
{
// If the loop propagator already exists, don't make the module again.
if (!(Environment::getInstance().hasModule(propName)))
{
MLoop::NoiseLoop::Par loopPar;
loopPar.q = resName;
loopPar.eta = srcName;
application.createModule<MLoop::NoiseLoop>(propName, loopPar);
}
}
/*******************************************************************************
* Contraction module creation.
******************************************************************************/
/*******************************************************************************
* Name: mesonContraction
* Purpose: Create meson contraction module and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* label - unique label to construct module name.
* mom - momentum to project (default is zero)
* gammas - gamma insertions at source and sink.
* Returns: None.
******************************************************************************/
inline void mesonContraction(Application &application, unsigned int npt,
std::string &q1, std::string &q2,
std::string &label,
std::string mom = ZERO_MOM,
std::string gammas = "<Gamma5 Gamma5>")
{
std::string modName = std::to_string(npt) + "pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::Meson::Par mesPar;
mesPar.output = std::to_string(npt) + "pt/" + label;
mesPar.q1 = q1;
mesPar.q2 = q2;
mesPar.mom = mom;
mesPar.gammas = gammas;
application.createModule<MContraction::Meson>(modName, mesPar);
}
}
/*******************************************************************************
* Name: gamma3ptContraction
* Purpose: Create gamma3pt contraction module and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* label - unique label to construct module name.
* gamma - gamma insertions between q2 and q3.
* Returns: None.
******************************************************************************/
inline void gamma3ptContraction(Application &application, unsigned int npt,
std::string &q1, std::string &q2,
std::string &q3, std::string &label,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
std::string modName = std::to_string(npt) + "pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::Gamma3pt::Par gamma3ptPar;
gamma3ptPar.output = std::to_string(npt) + "pt/" + label;
gamma3ptPar.q1 = q1;
gamma3ptPar.q2 = q2;
gamma3ptPar.q3 = q3;
gamma3ptPar.gamma = gamma;
application.createModule<MContraction::Gamma3pt>(modName, gamma3ptPar);
}
}
/*******************************************************************************
* Name: weakContraction[Eye,NonEye]
* Purpose: Create Weak Hamiltonian contraction module for Eye/NonEye topology
* and add to application module.
* Parameters: application - main application that stores modules.
* npt - specify n-point correlator (for labelling).
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* q4 - quark propagator 4.
* label - unique label to construct module name.
* Returns: None.
******************************************************************************/
#define HW_CONTRACTION(top) \
inline void weakContraction##top(Application &application, unsigned int npt,\
std::string &q1, std::string &q2, \
std::string &q3, std::string &q4, \
std::string &label)\
{\
std::string modName = std::to_string(npt) + "pt_" + label;\
if (!(Environment::getInstance().hasModule(modName)))\
{\
MContraction::WeakHamiltonian##top::Par weakPar;\
weakPar.output = std::to_string(npt) + "pt/" + label;\
weakPar.q1 = q1;\
weakPar.q2 = q2;\
weakPar.q3 = q3;\
weakPar.q4 = q4;\
application.createModule<MContraction::WeakHamiltonian##top>(modName, weakPar);\
}\
}
HW_CONTRACTION(Eye) // weakContractionEye
HW_CONTRACTION(NonEye) // weakContractionNonEye
/*******************************************************************************
* Name: disc0Contraction
* Purpose: Create contraction module for 4pt Weak Hamiltonian + current
* disconnected topology for neutral mesons and add to application
* module.
* Parameters: application - main application that stores modules.
* q1 - quark propagator 1.
* q2 - quark propagator 2.
* q3 - quark propagator 3.
* q4 - quark propagator 4.
* label - unique label to construct module name.
* Returns: None.
******************************************************************************/
inline void disc0Contraction(Application &application,
std::string &q1, std::string &q2,
std::string &q3, std::string &q4,
std::string &label)
{
std::string modName = "4pt_" + label;
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::WeakNeutral4ptDisc::Par disc0Par;
disc0Par.output = "4pt/" + label;
disc0Par.q1 = q1;
disc0Par.q2 = q2;
disc0Par.q3 = q3;
disc0Par.q4 = q4;
application.createModule<MContraction::WeakNeutral4ptDisc>(modName, disc0Par);
}
}
/*******************************************************************************
* Name: discLoopContraction
* Purpose: Create contraction module for disconnected loop and add to
* application module.
* Parameters: application - main application that stores modules.
* q_loop - loop quark propagator.
* modName - unique module name.
* gamma - gamma matrix to use in contraction.
* Returns: None.
******************************************************************************/
inline void discLoopContraction(Application &application,
std::string &q_loop, std::string &modName,
Gamma::Algebra gamma = Gamma::Algebra::Identity)
{
if (!(Environment::getInstance().hasModule(modName)))
{
MContraction::DiscLoop::Par discPar;
discPar.output = "disc/" + modName;
discPar.q_loop = q_loop;
discPar.gamma = gamma;
application.createModule<MContraction::DiscLoop>(modName, discPar);
}
}

22
tests/hadrons/Test_hadrons_meson_3pt.cc
View File

@ -30,14 +30,6 @@
using namespace Grid;
using namespace Hadrons;
static Gamma::Algebra gmu[4] =
{
Gamma::Algebra::GammaX,
Gamma::Algebra::GammaY,
Gamma::Algebra::GammaZ,
Gamma::Algebra::GammaT
};
int main(int argc, char *argv[])
{
// initialization //////////////////////////////////////////////////////////
@ -110,7 +102,7 @@ int main(int argc, char *argv[])
seqName.push_back(std::vector<std::string>(Nd));
for (unsigned int mu = 0; mu < Nd; ++mu)
{
seqPar.gamma = gmu[mu];
seqPar.gamma = 0x1 << mu;
seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
+ "_" + std::to_string(seqPar.tA) + "-"
+ qName[i];
@ -135,11 +127,11 @@ int main(int argc, char *argv[])
for (unsigned int i = 0; i < flavour.size(); ++i)
for (unsigned int j = i; j < flavour.size(); ++j)
{
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
mesPar.q1 = qName[i];
mesPar.q2 = qName[j];
mesPar.gammaSource = Gamma::Algebra::Gamma5;
mesPar.gammaSink = Gamma::Algebra::Gamma5;
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
mesPar.q1 = qName[i];
mesPar.q2 = qName[j];
mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_Z2_"
+ std::to_string(t)
+ "_"
@ -157,6 +149,8 @@ int main(int argc, char *argv[])
+ std::to_string(mu);
mesPar.q1 = qName[i];
mesPar.q2 = seqName[j][mu];
mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("3pt_Z2_"
+ std::to_string(t)
+ "_"

337
tests/hadrons/Test_hadrons_rarekaon.cc Normal file
View File

@ -0,0 +1,337 @@
/*******************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: tests/hadrons/Test_hadrons_rarekaon.cc
Copyright (C) 2017
Author: Andrew Lawson <andrew.lawson1991@gmail.com>
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.
*******************************************************************************/
#include "Test_hadrons.hpp"
using namespace Grid;
using namespace Hadrons;
enum quarks
{
light = 0,
strange = 1,
charm = 2
};
int main(int argc, char *argv[])
{
// parse command line //////////////////////////////////////////////////////
std::string configStem;
if (argc < 2)
{
std::cerr << "usage: " << argv[0] << " <configuration filestem> [Grid options]";
std::cerr << std::endl;
std::exit(EXIT_FAILURE);
}
configStem = argv[1];
// initialization //////////////////////////////////////////////////////////
Grid_init(&argc, &argv);
HadronsLogError.Active(GridLogError.isActive());
HadronsLogWarning.Active(GridLogWarning.isActive());
HadronsLogMessage.Active(GridLogMessage.isActive());
HadronsLogIterative.Active(GridLogIterative.isActive());
HadronsLogDebug.Active(GridLogDebug.isActive());
LOG(Message) << "Grid initialized" << std::endl;
// run setup ///////////////////////////////////////////////////////////////
Application application;
std::vector<double> mass = {.01, .04, .2};
std::vector<std::string> flavour = {"l", "s", "c"};
std::vector<std::string> solvers = {"CG_l", "CG_s", "CG_c"};
std::string kmom = "0. 0. 0. 0.";
std::string pmom = "1. 0. 0. 0.";
std::string qmom = "-1. 0. 0. 0.";
std::string mqmom = "1. 0. 0. 0.";
std::vector<unsigned int> tKs = {0};
unsigned int dt_pi = 16;
std::vector<unsigned int> tJs = {8};
unsigned int n_noise = 1;
unsigned int nt = 32;
bool do_disconnected(false);
// Global parameters.
Application::GlobalPar globalPar;
globalPar.trajCounter.start = 1500;
globalPar.trajCounter.end = 1520;
globalPar.trajCounter.step = 20;
globalPar.seed = "1 2 3 4";
globalPar.genetic.maxGen = 1000;
globalPar.genetic.maxCstGen = 200;
globalPar.genetic.popSize = 20;
globalPar.genetic.mutationRate = .1;
application.setPar(globalPar);
// gauge field
if (configStem == "None")
{
application.createModule<MGauge::Unit>("gauge");
}
else
{
MGauge::Load::Par gaugePar;
gaugePar.file = configStem;
application.createModule<MGauge::Load>("gauge", gaugePar);
}
for (unsigned int i = 0; i < flavour.size(); ++i)
{
// actions
MAction::DWF::Par actionPar;
actionPar.gauge = "gauge";
actionPar.Ls = 16;
actionPar.M5 = 1.8;
actionPar.mass = mass[i];
application.createModule<MAction::DWF>("DWF_" + flavour[i], actionPar);
// solvers
// RBPrecCG -> CG
MSolver::RBPrecCG::Par solverPar;
solverPar.action = "DWF_" + flavour[i];
solverPar.residual = 1.0e-8;
application.createModule<MSolver::RBPrecCG>(solvers[i],
solverPar);
}
// Create noise propagators for loops.
std::vector<std::string> noiseSrcs;
std::vector<std::vector<std::string>> noiseRes;
std::vector<std::vector<std::string>> noiseProps;
if (n_noise > 0)
{
MSource::Z2::Par noisePar;
noisePar.tA = 0;
noisePar.tB = nt - 1;
std::string loop_stem = "loop_";
noiseRes.resize(flavour.size());
noiseProps.resize(flavour.size());
for (unsigned int nn = 0; nn < n_noise; ++nn)
{
std::string eta = INIT_INDEX("noise", nn);
application.createModule<MSource::Z2>(eta, noisePar);
noiseSrcs.push_back(eta);
for (unsigned int f = 0; f < flavour.size(); ++f)
{
std::string loop_prop = INIT_INDEX(loop_stem + flavour[f], nn);
std::string loop_res = loop_prop + "_res";
makePropagator(application, loop_res, eta, solvers[f]);
makeLoop(application, loop_prop, eta, loop_res);
noiseRes[f].push_back(loop_res);
noiseProps[f].push_back(loop_prop);
}
}
}
// Translate rare kaon decay across specified timeslices.
for (unsigned int i = 0; i < tKs.size(); ++i)
{
// Zero-momentum wall source propagators for kaon and pion.
unsigned int tK = tKs[i];
unsigned int tpi = (tK + dt_pi) % nt;
std::string q_Kl_0 = INIT_INDEX("Q_l_0", tK);
std::string q_pil_0 = INIT_INDEX("Q_l_0", tpi);
MAKE_WALL_PROP(tK, q_Kl_0, solvers[light]);
MAKE_WALL_PROP(tpi, q_pil_0, solvers[light]);
// Wall sources for kaon and pion with momentum insertion. If either
// p or k are zero, or p = k, re-use the existing name to avoid
// duplicating a propagator.
std::string q_Ks_k = INIT_INDEX("Q_Ks_k", tK);
std::string q_Ks_p = INIT_INDEX((kmom == pmom) ? "Q_Ks_k" : "Q_Ks_p", tK);
std::string q_pil_k = INIT_INDEX((kmom == ZERO_MOM) ? "Q_l_0" : "Q_l_k", tpi);
std::string q_pil_p = INIT_INDEX((pmom == kmom) ? q_pil_k : ((pmom == ZERO_MOM) ? "Q_l_0" : "Q_l_p"), tpi);
MAKE_3MOM_WALL_PROP(tK, kmom, q_Ks_k, solvers[strange]);
MAKE_3MOM_WALL_PROP(tK, pmom, q_Ks_p, solvers[strange]);
MAKE_3MOM_WALL_PROP(tpi, kmom, q_pil_k, solvers[light]);
MAKE_3MOM_WALL_PROP(tpi, pmom, q_pil_p, solvers[light]);
/***********************************************************************
* CONTRACTIONS: pi and K 2pt contractions with mom = p, k.
**********************************************************************/
// Wall-Point
std::string PW_K_k = INIT_INDEX("PW_K_k", tK);
std::string PW_K_p = INIT_INDEX("PW_K_p", tK);
std::string PW_pi_k = INIT_INDEX("PW_pi_k", tpi);
std::string PW_pi_p = INIT_INDEX("PW_pi_p", tpi);
mesonContraction(application, 2, q_Kl_0, q_Ks_k, PW_K_k, kmom);
mesonContraction(application, 2, q_Kl_0, q_Ks_p, PW_K_p, pmom);
mesonContraction(application, 2, q_pil_k, q_pil_0, PW_pi_k, kmom);
mesonContraction(application, 2, q_pil_p, q_pil_0, PW_pi_p, pmom);
// Wall-Wall, to be done - requires modification of meson module.
/***********************************************************************
* CONTRACTIONS: 3pt Weak Hamiltonian, C & W (non-Eye type) classes.
**********************************************************************/
std::string HW_CW_k = LABEL_3PT("HW_CW_k", tK, tpi);
std::string HW_CW_p = LABEL_3PT("HW_CW_p", tK, tpi);
weakContractionNonEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, q_pil_0, HW_CW_k);
weakContractionNonEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, q_pil_0, HW_CW_p);
/***********************************************************************
* CONTRACTIONS: 3pt sd insertion.
**********************************************************************/
// Note: eventually will use wall sink smeared q_Kl_0 instead.
std::string sd_k = LABEL_3PT("sd_k", tK, tpi);
std::string sd_p = LABEL_3PT("sd_p", tK, tpi);
gamma3ptContraction(application, 3, q_Kl_0, q_Ks_k, q_pil_k, sd_k);
gamma3ptContraction(application, 3, q_Kl_0, q_Ks_p, q_pil_p, sd_p);
for (unsigned int nn = 0; nn < n_noise; ++nn)
{
/*******************************************************************
* CONTRACTIONS: 3pt Weak Hamiltonian, S and E (Eye type) classes.
******************************************************************/
// Note: eventually will use wall sink smeared q_Kl_0 instead.
for (unsigned int f = 0; f < flavour.size(); ++f)
{
if ((f != strange) || do_disconnected)
{
std::string HW_SE_k = LABEL_3PT("HW_SE_k_" + flavour[f], tK, tpi);
std::string HW_SE_p = LABEL_3PT("HW_SE_p_" + flavour[f], tK, tpi);
std::string loop_q = noiseProps[f][nn];
weakContractionEye(application, 3, q_Kl_0, q_Ks_k, q_pil_k, loop_q, HW_CW_k);
weakContractionEye(application, 3, q_Kl_0, q_Ks_p, q_pil_p, loop_q, HW_CW_p);
}
}
}
// Perform separate contractions for each t_J position.
for (unsigned int j = 0; j < tJs.size(); ++j)
{
// Sequential sources for current insertions. Local for now,
// gamma_0 only.
unsigned int tJ = (tJs[j] + tK) % nt;
MSource::SeqGamma::Par seqPar;
std::string q_KlCl_q = LABEL_3PT("Q_KlCl_q", tK, tJ);
std::string q_KsCs_mq = LABEL_3PT("Q_KsCs_mq", tK, tJ);
std::string q_pilCl_q = LABEL_3PT("Q_pilCl_q", tpi, tJ);
std::string q_pilCl_mq = LABEL_3PT("Q_pilCl_mq", tpi, tJ);
MAKE_SEQUENTIAL_PROP(tJ, q_Kl_0, qmom, q_KlCl_q, solvers[light]);
MAKE_SEQUENTIAL_PROP(tJ, q_Ks_k, mqmom, q_KsCs_mq, solvers[strange]);
MAKE_SEQUENTIAL_PROP(tJ, q_pil_p, qmom, q_pilCl_q, solvers[light]);
MAKE_SEQUENTIAL_PROP(tJ, q_pil_0, mqmom, q_pilCl_mq, solvers[light]);
/*******************************************************************
* CONTRACTIONS: pi and K 3pt contractions with current insertion.
******************************************************************/
// Wall-Point
std::string C_PW_Kl = LABEL_3PT("C_PW_Kl", tK, tJ);
std::string C_PW_Ksb = LABEL_3PT("C_PW_Ksb", tK, tJ);
std::string C_PW_pilb = LABEL_3PT("C_PW_pilb", tK, tJ);
std::string C_PW_pil = LABEL_3PT("C_PW_pil", tK, tJ);
mesonContraction(application, 3, q_KlCl_q, q_Ks_k, C_PW_Kl, pmom);
mesonContraction(application, 3, q_Kl_0, q_KsCs_mq, C_PW_Ksb, pmom);
mesonContraction(application, 3, q_pil_0, q_pilCl_q, C_PW_pilb, kmom);
mesonContraction(application, 3, q_pilCl_mq, q_pil_p, C_PW_pil, kmom);
// Wall-Wall, to be done.
/*******************************************************************
* CONTRACTIONS: 4pt contractions, C & W classes.
******************************************************************/
std::string CW_Kl = LABEL_4PT("CW_Kl", tK, tJ, tpi);
std::string CW_Ksb = LABEL_4PT("CW_Ksb", tK, tJ, tpi);
std::string CW_pilb = LABEL_4PT("CW_pilb", tK, tJ, tpi);
std::string CW_pil = LABEL_4PT("CW_pil", tK, tJ, tpi);
weakContractionNonEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, q_pil_0, CW_Kl);
weakContractionNonEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, q_pil_0, CW_Ksb);
weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, q_pil_0, CW_pilb);
weakContractionNonEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, q_pilCl_mq, CW_pil);
/*******************************************************************
* CONTRACTIONS: 4pt contractions, sd insertions.
******************************************************************/
// Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
std::string sd_Kl = LABEL_4PT("sd_Kl", tK, tJ, tpi);
std::string sd_Ksb = LABEL_4PT("sd_Ksb", tK, tJ, tpi);
std::string sd_pilb = LABEL_4PT("sd_pilb", tK, tJ, tpi);
gamma3ptContraction(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, sd_Kl);
gamma3ptContraction(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, sd_Ksb);
gamma3ptContraction(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, sd_pilb);
// Sequential sources for each noise propagator.
for (unsigned int nn = 0; nn < n_noise; ++nn)
{
std::string loop_stem = "loop_";
// Contraction required for each quark flavour - alternatively
// drop the strange loop if not performing disconnected
// contractions or neglecting H_W operators Q_3 -> Q_10.
for (unsigned int f = 0; f < flavour.size(); ++f)
{
if ((f != strange) || do_disconnected)
{
std::string eta = noiseSrcs[nn];
std::string loop_q = noiseProps[f][nn];
std::string loop_qCq = LABEL_3PT(loop_stem + flavour[f], tJ, nn);
std::string loop_qCq_res = loop_qCq + "_res";
MAKE_SEQUENTIAL_PROP(tJ, noiseRes[f][nn], qmom,
loop_qCq_res, solvers[f]);
makeLoop(application, loop_qCq, eta, loop_qCq_res);
/*******************************************************
* CONTRACTIONS: 4pt contractions, S & E classes.
******************************************************/
// Note: eventually will use wall sink smeared q_Kl_0/q_KlCl_q instead.
std::string SE_Kl = LABEL_4PT_NOISE("SE_Kl", tK, tJ, tpi, nn);
std::string SE_Ksb = LABEL_4PT_NOISE("SE_Ksb", tK, tJ, tpi, nn);
std::string SE_pilb = LABEL_4PT_NOISE("SE_pilb", tK, tJ, tpi, nn);
std::string SE_loop = LABEL_4PT_NOISE("SE_loop", tK, tJ, tpi, nn);
weakContractionEye(application, 4, q_KlCl_q, q_Ks_k, q_pil_p, loop_q, SE_Kl);
weakContractionEye(application, 4, q_Kl_0, q_KsCs_mq, q_pil_p, loop_q, SE_Ksb);
weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, SE_pilb);
weakContractionEye(application, 4, q_Kl_0, q_Ks_k, q_pil_p, loop_qCq, SE_loop);
/*******************************************************
* CONTRACTIONS: 4pt contractions, pi0 disconnected
* loop.
******************************************************/
std::string disc0 = LABEL_4PT_NOISE("disc0", tK, tJ, tpi, nn);
disc0Contraction(application, q_Kl_0, q_Ks_k, q_pilCl_q, loop_q, disc0);
/*******************************************************
* CONTRACTIONS: Disconnected loop.
******************************************************/
std::string discLoop = "disc_" + loop_qCq;
discLoopContraction(application, loop_qCq, discLoop);
}
}
}
}
}
// execution
std::string par_file_name = "rarekaon_000_100_tK0_tpi16_tJ8_noloop_mc0.2.xml";
application.saveParameterFile(par_file_name);
application.run();
// epilogue
LOG(Message) << "Grid is finalizing now" << std::endl;
Grid_finalize();
return EXIT_SUCCESS;
}

4
tests/hadrons/Test_hadrons_spectrum.cc
View File

@ -96,12 +96,16 @@ int main(int argc, char *argv[])
mesPar.output = "mesons/pt_" + flavour[i] + flavour[j];
mesPar.q1 = "Qpt_" + flavour[i];
mesPar.q2 = "Qpt_" + flavour[j];
mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_pt_"
+ flavour[i] + flavour[j],
mesPar);
mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
mesPar.q1 = "QZ2_" + flavour[i];
mesPar.q2 = "QZ2_" + flavour[j];
mesPar.gammas = "all";
mesPar.mom = "0. 0. 0. 0.";
application.createModule<MContraction::Meson>("meson_Z2_"
+ flavour[i] + flavour[j],
mesPar);

285
tests/hmc/Test_hmc_EOMobiusRatio.cc Normal file
View File

@ -0,0 +1,285 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_hmc_EODWFRatio.cc
Copyright (C) 2015-2016
Author: Peter Boyle <pabobyle@ph.ed.ac.uk>
Author: Guido Cossu <guido.cossu@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>
namespace Grid{
struct FermionParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(FermionParameters,
int, Ls,
double, mass,
double, M5,
double, b,
double, c,
double, StoppingCondition,
int, MaxCGIterations,
bool, ApplySmearing);
//template <class ReaderClass >
//FermionParameters(Reader<ReaderClass>& Reader){
// read(Reader, "Mobius", *this);
//}
};
struct MobiusHMCParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(MobiusHMCParameters,
double, gauge_beta,
FermionParameters, Mobius)
template <class ReaderClass >
MobiusHMCParameters(Reader<ReaderClass>& Reader){
read(Reader, "Action", *this);
}
};
struct SmearingParameters: Serializable {
GRID_SERIALIZABLE_CLASS_MEMBERS(SmearingParameters,
double, rho,
Integer, Nsmear)
template <class ReaderClass >
SmearingParameters(Reader<ReaderClass>& Reader){
read(Reader, "StoutSmearing", *this);
}
};
}
int main(int argc, char **argv) {
using namespace Grid;
using namespace Grid::QCD;
Grid_init(&argc, &argv);
int threads = GridThread::GetThreads();
// here make a routine to print all the relevant information on the run
std::cout << GridLogMessage << "Grid is setup to use " << threads << " threads" << std::endl;
// Typedefs to simplify notation
typedef GenericHMCRunner<MinimumNorm2> HMCWrapper; // Uses the default minimum norm
typedef WilsonImplR FermionImplPolicy;
typedef MobiusFermionR FermionAction;
typedef typename FermionAction::FermionField FermionField;
// Serialiser
//typedef Grid::XmlReader Serialiser;
typedef Grid::JSONReader Serialiser;
//::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
HMCWrapper TheHMC;
TheHMC.ReadCommandLine(argc, argv); // these can be parameters from file
// Reader, file should come from command line
if (TheHMC.ParameterFile.empty()){
std::cout << "Input file not specified."
<< "Use --ParameterFile option in the command line.\nAborting"
<< std::endl;
exit(1);
}
Serialiser Reader(TheHMC.ParameterFile);
MobiusHMCParameters MyParams(Reader);
// Apply smearing to the fermionic action
bool ApplySmearing = MyParams.Mobius.ApplySmearing;
// Grid from the command line
TheHMC.Resources.AddFourDimGrid("gauge");
// Possibile to create the module by hand
// hardcoding parameters or using a Reader
// Checkpointer definition (Name: Checkpointer)
CheckpointerParameters CPparams(Reader);
// Commenting out since we are using the reader
/*
CPparams.config_prefix = "ckpoint_EODWF_lat";
CPparams.rng_prefix = "ckpoint_EODWF_rng";
CPparams.saveInterval = 5;
CPparams.format = "IEEE64BIG";
*/
TheHMC.Resources.LoadBinaryCheckpointer(CPparams);
// RNG definition (Name: RandomNumberGenerator)
RNGModuleParameters RNGpar(Reader);
// Commenting out since we are using the reader
/*
RNGpar.serial_seeds = "1 2 3 4 5";
RNGpar.parallel_seeds = "6 7 8 9 10";
*/
TheHMC.Resources.SetRNGSeeds(RNGpar);
// Construct observables
// here there is too much indirection
typedef PlaquetteMod<HMCWrapper::ImplPolicy> PlaqObs;
TheHMC.Resources.AddObservable<PlaqObs>();
//////////////////////////////////////////////
/////////////////////////////////////////////////////////////
// Collect actions, here use more encapsulation
// need wrappers of the fermionic classes
// that have a complex construction
// standard
//RealD beta = 5.6 ;
WilsonGaugeActionR Waction(MyParams.gauge_beta);
//const int Ls = 8;
const int Ls = MyParams.Mobius.Ls;
auto GridPtr = TheHMC.Resources.GetCartesian();
auto GridRBPtr = TheHMC.Resources.GetRBCartesian();
auto FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,GridPtr);
auto FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,GridPtr);
// temporarily need a gauge field
LatticeGaugeField U(GridPtr);
Real mass = MyParams.Mobius.mass; //0.04;
Real pv = 1.0;
RealD M5 = MyParams.Mobius.M5; //1.5;
// Note: IroIro and Grid notation for b and c differ
RealD b = MyParams.Mobius.b; // 3./2.;
RealD c = MyParams.Mobius.c; // 1./2.;
// These lines are unecessary if BC are all periodic
std::vector<Complex> boundary = {1,1,1,-1};
FermionAction::ImplParams Params(boundary);
FermionAction DenOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,mass,M5,b,c, Params);
FermionAction NumOp(U,*FGrid,*FrbGrid,*GridPtr,*GridRBPtr,pv, M5,b,c, Params);
//double StoppingCondition = 1e-8;
//double MaxCGIterations = 10000;
ConjugateGradient<FermionField> CG(MyParams.Mobius.StoppingCondition,MyParams.Mobius.MaxCGIterations);
TwoFlavourEvenOddRatioPseudoFermionAction<FermionImplPolicy> Nf2(NumOp, DenOp,CG,CG);
// Set smearing (true/false), default: false
Nf2.is_smeared = ApplySmearing;
// Collect actions
ActionLevel<HMCWrapper::Field> Level1(1);
Level1.push_back(&Nf2);
ActionLevel<HMCWrapper::Field> Level2(4);
Level2.push_back(&Waction);
TheHMC.TheAction.push_back(Level1);
TheHMC.TheAction.push_back(Level2);
/////////////////////////////////////////////////////////////
// HMC parameters are serialisable
TheHMC.Parameters.initialize(Reader);
/*
TheHMC.Parameters.MD.MDsteps = 20;
TheHMC.Parameters.MD.trajL = 1.0;
*/
// Reset performance counters
NumOp.ZeroCounters();
DenOp.ZeroCounters();
if (ApplySmearing){
SmearingParameters SmPar(Reader);
//double rho = 0.1; // smearing parameter
//int Nsmear = 3; // number of smearing levels
Smear_Stout<HMCWrapper::ImplPolicy> Stout(SmPar.rho);
SmearedConfiguration<HMCWrapper::ImplPolicy> SmearingPolicy(GridPtr, SmPar.Nsmear, Stout);
TheHMC.Run(SmearingPolicy); // for smearing
} else {
TheHMC.Run(); // no smearing
}
std::cout << GridLogMessage << "Numerator report, Pauli-Villars term : " << std::endl;
NumOp.Report();
std::cout << GridLogMessage << "Denominator report, Dw(m) term (includes CG) : " << std::endl;
DenOp.Report();
Grid_finalize();
} // main
/* Examples for input files
JSON
{
"Checkpointer": {
"config_prefix": "ckpoint_json_lat",
"rng_prefix": "ckpoint_json_rng",
"saveInterval": 1,
"format": "IEEE64BIG"
},
"RandomNumberGenerator": {
"serial_seeds": "1 2 3 4 6",
"parallel_seeds": "6 7 8 9 11"
},
"Action":{
"gauge_beta": 5.6,
"Mobius": {
"Ls" : 10,
"mass": 0.01,
"M5" : 1.0,
"b" : 1.5,
"c" : 0.5,
"StoppingCondition": 1e-8,
"MaxCGIterations": 10000,
"ApplySmearing": true
}
},
"HMC":{
"StartTrajectory": 0,
"Trajectories": 100,
"MetropolisTest": true,
"NoMetropolisUntil": 10,
"StartingType": "HotStart",
"MD":{
"name": "MinimumNorm2",
"MDsteps": 15,
"trajL": 2.0
}
},
"StoutSmearing":{
"rho": 0.1,
"Nsmear": 3
}
}
XML example not provided yet
*/

44
tests/qdpxx/Test_qdpxx_loops_staples.cc
View File

@ -317,6 +317,7 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
Grid::QCD::LatticeComplex rect(UGrid);
Grid::QCD::TComplex trect;
Grid::QCD::Complex crect;
Grid::RealD rrect;
Grid::RealD vol = UGrid->gSites();
for(int mu=0;mu<Grid::QCD::Nd;mu++){
for(int nu=0;nu<Grid::QCD::Nd;nu++){
@ -325,7 +326,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
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;
rrect = real(crect);
std::cout<< "mu/nu = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
Grid::GridStopWatch Peter;
Grid::GridStopWatch Azusa;
@ -355,7 +357,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect=real(crect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __
@ -370,7 +373,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect=real(crect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
@ -384,7 +388,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __ ___
@ -399,7 +404,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// --
@ -423,7 +429,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
SumTrStap += TrStap;
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
@ -441,11 +448,13 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
TrStap = Grid::trace (U[mu]*Stap);
trect = Grid::sum(TrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
trect = Grid::sum(SumTrStap);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/2.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline trace 2x1+1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/2.0/3.0<<std::endl;
}
Peter.Stop();
Azusa.Start();
@ -489,7 +498,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __ __
// |__ |
@ -501,7 +511,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __
// |__ __ |
@ -513,7 +524,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// __
@ -526,7 +538,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 2x1 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// 1(mu) x 2 (nu) ** this was ok
@ -542,7 +555,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
// 1(mu) x 2 (nu) ** this was ok
//
@ -570,8 +584,8 @@ double calc_grid_r_dir(Grid::QCD::LatticeGaugeField & Umu)
RectPlaq_d = Grid::trace(U[mu]*ds_U);
trect = Grid::sum(RectPlaq_d);
crect = Grid::TensorRemove(trect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<crect/vol/1.0/3.0<<std::endl;
rrect = real(crect);
std::cout<< "mu/nu inline AZUSA trace 1x2 code = "<<mu<<"/"<<nu<<" ; rect = "<<rrect/vol/1.0/3.0<<std::endl;
}
Azusa.Stop();

2
tests/solver/Test_dwf_cg_prec.cc
View File

@ -89,7 +89,7 @@ int main(int argc, char** argv) {
GridStopWatch CGTimer;
SchurDiagMooeeOperator<DomainWallFermionR, LatticeFermion> HermOpEO(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8, 10000, 0);// switch off the assert
ConjugateGradient<LatticeFermion> CG(1.0e-5, 10000, 0);// switch off the assert
CGTimer.Start();
CG(HermOpEO, src_o, result_o);

2
tests/solver/Test_dwf_cg_unprec.cc
View File

@ -73,7 +73,7 @@ int main (int argc, char ** argv)
DomainWallFermionR Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
MdagMLinearOperator<DomainWallFermionR,LatticeFermion> HermOp(Ddwf);
ConjugateGradient<LatticeFermion> CG(1.0e-8,10000);
ConjugateGradient<LatticeFermion> CG(1.0e-6,10000);
CG(HermOp,src,result);
Grid_finalize();

119
tests/solver/Test_staggered_block_cg_unprec.cc Normal file
View File

@ -0,0 +1,119 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_cg_unprec.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
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>
#include <Grid/algorithms/iterative/BlockConjugateGradient.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)
{
typedef typename ImprovedStaggeredFermion5DR::FermionField FermionField;
typedef typename ImprovedStaggeredFermion5DR::ComplexField ComplexField;
typename ImprovedStaggeredFermion5DR::ImplParams params;
const int Ls=4;
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 * 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> seeds({1,2,3,4});
GridParallelRNG pRNG(UGrid ); pRNG.SeedFixedIntegers(seeds);
GridParallelRNG pRNG5(FGrid); pRNG5.SeedFixedIntegers(seeds);
FermionField src(FGrid); random(pRNG5,src);
FermionField result(FGrid); result=zero;
RealD nrm = norm2(src);
LatticeGaugeField Umu(UGrid); SU3::HotConfiguration(pRNG,Umu);
RealD mass=0.01;
ImprovedStaggeredFermion5DR Ds(Umu,Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass);
MdagMLinearOperator<ImprovedStaggeredFermion5DR,FermionField> HermOp(Ds);
ConjugateGradient<FermionField> CG(1.0e-8,10000);
BlockConjugateGradient<FermionField> BCG(1.0e-8,10000);
MultiRHSConjugateGradient<FermionField> mCG(1.0e-8,10000);
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << " Calling 4d CG "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
ImprovedStaggeredFermionR Ds4d(Umu,Umu,*UGrid,*UrbGrid,mass);
MdagMLinearOperator<ImprovedStaggeredFermionR,FermionField> HermOp4d(Ds4d);
FermionField src4d(UGrid); random(pRNG,src4d);
FermionField result4d(UGrid); result4d=zero;
CG(HermOp4d,src4d,result4d);
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << " Calling 5d CG for "<<Ls <<" right hand sides" <<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
result=zero;
CG(HermOp,src,result);
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << " Calling multiRHS CG for "<<Ls <<" right hand sides" <<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
result=zero;
mCG(HermOp,src,result);
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
std::cout << GridLogMessage << " Calling Block CG for "<<Ls <<" right hand sides" <<std::endl;
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
result=zero;
BCG(HermOp,src,result);
std::cout << GridLogMessage << "************************************************************************ "<<std::endl;
Grid_finalize();
}

83
tests/solver/Test_staggered_cg_unprec.cc Normal file
View File

@ -0,0 +1,83 @@
/*************************************************************************************
Grid physics library, www.github.com/paboyle/Grid
Source file: ./tests/Test_wilson_cg_unprec.cc
Copyright (C) 2015
Author: Azusa Yamaguchi <ayamaguc@staffmail.ed.ac.uk>
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>
#include <Grid/algorithms/iterative/BlockConjugateGradient.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)
{
typedef typename ImprovedStaggeredFermionR::FermionField FermionField;
typedef typename ImprovedStaggeredFermionR::ComplexField ComplexField;
typename ImprovedStaggeredFermionR::ImplParams params;
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);
std::vector<int> seeds({1,2,3,4});
GridParallelRNG pRNG(&Grid); pRNG.SeedFixedIntegers(seeds);
FermionField src(&Grid); random(pRNG,src);
RealD nrm = norm2(src);
FermionField result(&Grid); result=zero;
LatticeGaugeField Umu(&Grid); SU3::HotConfiguration(pRNG,Umu);
double volume=1;
for(int mu=0;mu<Nd;mu++){
volume=volume*latt_size[mu];
}
RealD mass=0.1;
ImprovedStaggeredFermionR Ds(Umu,Umu,Grid,RBGrid,mass);
MdagMLinearOperator<ImprovedStaggeredFermionR,FermionField> HermOp(Ds);
ConjugateGradient<FermionField> CG(1.0e-6,10000);
CG(HermOp,src,result);
Grid_finalize();
}