mirror of
https://github.com/paboyle/Grid.git
synced 2025-06-14 13:57:07 +01:00
Merge branch 'master' of https://github.com/paboyle/Grid
This commit is contained in:
Azusa Yamaguchi
0
tests/InvSqrt.gnu
Normal file
0
tests/InvSqrt.gnu
Normal file
@ -1,5 +1,5 @@
|
||||
|
||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_even_odd Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_even_odd Test_gamma Test_main Test_multishift_sqrt Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_even_odd
|
||||
bin_PROGRAMS = Test_GaugeAction Test_cayley_cg Test_cayley_coarsen_support Test_cayley_even_odd Test_cayley_ldop_cg Test_cayley_ldop_cr Test_cf_coarsen_support Test_cf_cr_unprec Test_contfrac_cg Test_contfrac_even_odd Test_cshift Test_cshift_red_black Test_dwf_cg_prec Test_dwf_cg_schur Test_dwf_cg_unprec Test_dwf_cr_unprec Test_dwf_even_odd Test_gamma Test_main Test_multishift_sqrt Test_nersc_io Test_remez Test_rng Test_rng_fixed Test_simd Test_stencil Test_wilson_cg_prec Test_wilson_cg_schur Test_wilson_cg_unprec Test_wilson_cr_unprec Test_wilson_even_odd
|
||||
|
||||
|
||||
Test_GaugeAction_SOURCES=Test_GaugeAction.cc
|
||||
@ -10,10 +10,30 @@ Test_cayley_cg_SOURCES=Test_cayley_cg.cc
|
||||
Test_cayley_cg_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_coarsen_support_SOURCES=Test_cayley_coarsen_support.cc
|
||||
Test_cayley_coarsen_support_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_even_odd_SOURCES=Test_cayley_even_odd.cc
|
||||
Test_cayley_even_odd_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_ldop_cg_SOURCES=Test_cayley_ldop_cg.cc
|
||||
Test_cayley_ldop_cg_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cayley_ldop_cr_SOURCES=Test_cayley_ldop_cr.cc
|
||||
Test_cayley_ldop_cr_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cf_coarsen_support_SOURCES=Test_cf_coarsen_support.cc
|
||||
Test_cf_coarsen_support_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_cf_cr_unprec_SOURCES=Test_cf_cr_unprec.cc
|
||||
Test_cf_cr_unprec_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_contfrac_cg_SOURCES=Test_contfrac_cg.cc
|
||||
Test_contfrac_cg_LDADD=-lGrid
|
||||
|
||||
@ -42,6 +62,10 @@ Test_dwf_cg_unprec_SOURCES=Test_dwf_cg_unprec.cc
|
||||
Test_dwf_cg_unprec_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_dwf_cr_unprec_SOURCES=Test_dwf_cr_unprec.cc
|
||||
Test_dwf_cr_unprec_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_dwf_even_odd_SOURCES=Test_dwf_even_odd.cc
|
||||
Test_dwf_even_odd_LDADD=-lGrid
|
||||
|
||||
@ -94,6 +118,10 @@ Test_wilson_cg_unprec_SOURCES=Test_wilson_cg_unprec.cc
|
||||
Test_wilson_cg_unprec_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_wilson_cr_unprec_SOURCES=Test_wilson_cr_unprec.cc
|
||||
Test_wilson_cr_unprec_LDADD=-lGrid
|
||||
|
||||
|
||||
Test_wilson_even_odd_SOURCES=Test_wilson_even_odd.cc
|
||||
Test_wilson_even_odd_LDADD=-lGrid
|
||||
|
||||
|
||||
2
tests/Sqrt.gnu
Normal file
2
tests/Sqrt.gnu
Normal file
@ -0,0 +1,2 @@
|
||||
f(x) = 6.81384+(-2.34645e-06/(x+0.000228091))+(-1.51593e-05/(x+0.00112084))+(-6.89254e-05/(x+0.003496))+(-0.000288983/(x+0.00954309))+(-0.00119277/(x+0.024928))+(-0.0050183/(x+0.0646627))+(-0.0226449/(x+0.171576))+(-0.123767/(x+0.491792))+(-1.1705/(x+1.78667))+(-102.992/(x+18.4866));
|
||||
f(x) = 0.14676+(0.00952992/(x+5.40933e-05))+(0.0115952/(x+0.000559699))+(0.0161824/(x+0.00203338))+(0.0243252/(x+0.00582831))+(0.0379533/(x+0.0154649))+(0.060699/(x+0.0401156))+(0.100345/(x+0.104788))+(0.178335/(x+0.286042))+(0.381586/(x+0.892189))+(1.42625/(x+4.38422));
|
||||
@ -115,7 +115,7 @@ int main (int argc, char ** argv)
|
||||
Complex l = TensorRemove(Tl);
|
||||
std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
|
||||
|
||||
sumBlocks(cPlaq,Plaq);
|
||||
blockSum(cPlaq,Plaq);
|
||||
TComplex TcP = sum(cPlaq);
|
||||
Complex ll= TensorRemove(TcP);
|
||||
std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
|
||||
|
||||
174
tests/Test_cayley_coarsen_support.cc
Normal file
174
tests/Test_cayley_coarsen_support.cc
Normal file
@ -0,0 +1,174 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
|
||||
|
||||
#if 0
|
||||
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||
Umu=zero;
|
||||
Complex cone(1.0,0.0);
|
||||
for(int nn=0;nn<Nd;nn++){
|
||||
if(1) {
|
||||
if (nn>2) { U[nn]=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl; }
|
||||
else { U[nn]=cone; std::cout << "unit gauge field in dir "<<nn<<std::endl; }
|
||||
}
|
||||
pokeIndex<LorentzIndex>(Umu,U[nn],nn);
|
||||
}
|
||||
#endif
|
||||
|
||||
RealD mass=0.5;
|
||||
RealD M5=1.8;
|
||||
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
|
||||
HermIndefOp.Op(src,ref);
|
||||
HermIndefOp.OpDiag(src,result);
|
||||
|
||||
for(int d=0;d<4;d++){
|
||||
HermIndefOp.OpDir(src,tmp,d+1,+1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
HermIndefOp.OpDir(src,tmp,d+1,-1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
}
|
||||
err = result-ref;
|
||||
std::cout<<"Error "<<norm2(err)<<std::endl;
|
||||
|
||||
const int nbasis = 2;
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
LatticeFermion prom(FGrid);
|
||||
|
||||
for(int b=0;b<nbasis;b++){
|
||||
random(RNG5,subspace[b]);
|
||||
}
|
||||
std::cout << "Computed randoms"<< std::endl;
|
||||
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
LittleDiracOperator LittleDiracOp(*Coarse5d);
|
||||
|
||||
LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
CoarseVector c_proj(Coarse5d);
|
||||
|
||||
// TODO
|
||||
// -- promote from subspace, check we get the vector we wanted
|
||||
// -- apply ldop; check we get the same as inner product of M times big vec
|
||||
// -- pick blocks one by one. Evaluate matrix elements.
|
||||
Complex one(1.0);
|
||||
c_src = one; // 1 in every element for vector 1.
|
||||
|
||||
blockPromote(c_src,err,subspace);
|
||||
|
||||
|
||||
prom=zero;
|
||||
for(int b=0;b<nbasis;b++){
|
||||
prom=prom+subspace[b];
|
||||
}
|
||||
err=err-prom;
|
||||
std::cout<<"Promoted back from subspace err "<<norm2(err)<<std::endl;
|
||||
|
||||
HermIndefOp.HermOp(prom,tmp);
|
||||
blockProject(c_proj,tmp,subspace);
|
||||
|
||||
LittleDiracOp.M(c_src,c_res);
|
||||
|
||||
c_proj = c_proj - c_res;
|
||||
std::cout<<"Representation of ldop within subspace "<<norm2(c_proj)<<std::endl;
|
||||
|
||||
std::cout << "Multiplying by LittleDiracOp "<< std::endl;
|
||||
LittleDiracOp.M(c_src,c_res);
|
||||
|
||||
std::cout<<"Testing hermiticity explicitly by inspecting matrix elements"<<std::endl;
|
||||
LittleDiracOp.AssertHermitian();
|
||||
|
||||
std::cout << "Testing Hermiticity stochastically "<< std::endl;
|
||||
CoarseVector phi(Coarse5d);
|
||||
CoarseVector chi(Coarse5d);
|
||||
CoarseVector Aphi(Coarse5d);
|
||||
CoarseVector Achi(Coarse5d);
|
||||
|
||||
random(CRNG,phi);
|
||||
random(CRNG,chi);
|
||||
|
||||
|
||||
std::cout<<"Made randoms"<<std::endl;
|
||||
|
||||
LittleDiracOp.M(phi,Aphi);
|
||||
LittleDiracOp.Mdag(chi,Achi);
|
||||
|
||||
ComplexD pAc = innerProduct(chi,Aphi);
|
||||
ComplexD cAp = innerProduct(phi,Achi);
|
||||
ComplexD cAc = innerProduct(chi,Achi);
|
||||
ComplexD pAp = innerProduct(phi,Aphi);
|
||||
|
||||
std::cout<< "pAc "<<pAc<<" cAp "<< cAp<< " diff "<<pAc-adj(cAp)<<std::endl;
|
||||
|
||||
std::cout<< "pAp "<<pAp<<" cAc "<< cAc<<"Should be real"<< std::endl;
|
||||
|
||||
std::cout<<"Testing linearity"<<std::endl;
|
||||
CoarseVector PhiPlusChi(Coarse5d);
|
||||
CoarseVector APhiPlusChi(Coarse5d);
|
||||
CoarseVector linerr(Coarse5d);
|
||||
PhiPlusChi = phi+chi;
|
||||
LittleDiracOp.M(PhiPlusChi,APhiPlusChi);
|
||||
|
||||
linerr= APhiPlusChi-Aphi;
|
||||
linerr= linerr-Achi;
|
||||
std::cout<<"**Diff "<<norm2(linerr)<<std::endl;
|
||||
|
||||
|
||||
std::cout << "Done "<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
136
tests/Test_cayley_ldop_cg.cc
Normal file
136
tests/Test_cayley_ldop_cg.cc
Normal file
@ -0,0 +1,136 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
double lowpass(double x)
|
||||
{
|
||||
return pow(x*x+1.0,-2);
|
||||
}
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
|
||||
ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
|
||||
filter.csv(csv);
|
||||
csv.close();
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); gaussian(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu);
|
||||
|
||||
|
||||
//random(RNG4,Umu);
|
||||
NerscField header;
|
||||
std::string file("./ckpoint_lat.400");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
#if 0
|
||||
LatticeColourMatrix U(UGrid);
|
||||
U=zero;
|
||||
for(int nn=0;nn<Nd;nn++){
|
||||
if (nn>2) {
|
||||
pokeIndex<LorentzIndex>(Umu,U,nn);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.0;
|
||||
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
|
||||
const int nbasis = 8;
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
LatticeFermion noise(FGrid);
|
||||
LatticeFermion ms(FGrid);
|
||||
for(int b=0;b<nbasis;b++){
|
||||
Gamma g5(Gamma::Gamma5);
|
||||
gaussian(RNG5,noise);
|
||||
RealD scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
HermIndefOp.HermOp(noise,ms); std::cout << "Noise "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
|
||||
|
||||
// filter(HermIndefOp,noise,subspace[b]);
|
||||
// inverse iteration
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-5,10000);
|
||||
|
||||
for(int i=0;i<4;i++){
|
||||
|
||||
CG(HermDefOp,noise,subspace[b]);
|
||||
noise = subspace[b];
|
||||
|
||||
scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
HermDefOp.HermOp(noise,ms); std::cout << "filt "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
}
|
||||
|
||||
subspace[b] = noise;
|
||||
HermIndefOp.HermOp(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
|
||||
}
|
||||
std::cout << "Computed randoms"<< std::endl;
|
||||
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
LittleDiracOperator LittleDiracOp(*Coarse5d);
|
||||
LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
std::cout << "Solving CG on coarse space "<< std::endl;
|
||||
MdagMLinearOperator<LittleDiracOperator,CoarseVector> PosdefLdop(LittleDiracOp);
|
||||
ConjugateGradient<CoarseVector> CG(1.0e-8,10000);
|
||||
CG(PosdefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "Done "<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
138
tests/Test_cayley_ldop_cr.cc
Normal file
138
tests/Test_cayley_ldop_cr.cc
Normal file
@ -0,0 +1,138 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
double lowpass(double x)
|
||||
{
|
||||
return pow(x*x+1.0,-2);
|
||||
}
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
Chebyshev<LatticeFermion> filter(-150.0, 150.0,16, lowpass);
|
||||
ofstream csv(std::string("filter.dat"),std::ios::out|std::ios::trunc);
|
||||
filter.csv(csv);
|
||||
csv.close();
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
// Construct a coarsened grid
|
||||
std::vector<int> clatt = GridDefaultLatt();
|
||||
for(int d=0;d<clatt.size();d++){
|
||||
clatt[d] = clatt[d]/2;
|
||||
}
|
||||
GridCartesian *Coarse4d = SpaceTimeGrid::makeFourDimGrid(clatt, GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());;
|
||||
GridCartesian *Coarse5d = SpaceTimeGrid::makeFiveDimGrid(1,Coarse4d);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
std::vector<int> cseeds({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
GridParallelRNG CRNG(Coarse5d);CRNG.SeedFixedIntegers(cseeds);
|
||||
|
||||
LatticeFermion src(FGrid); gaussian(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid); gaussian(RNG4,Umu);
|
||||
|
||||
|
||||
//random(RNG4,Umu);
|
||||
NerscField header;
|
||||
std::string file("./ckpoint_lat.400");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
#if 0
|
||||
LatticeColourMatrix U(UGrid);
|
||||
Complex cone(1.0,0.0);
|
||||
for(int nn=0;nn<Nd;nn++){
|
||||
if (nn==3) {
|
||||
U=zero; std::cout << "zeroing gauge field in dir "<<nn<<std::endl;
|
||||
// else { U[nn]= cone;std::cout << "unit gauge field in dir "<<nn<<std::endl; }
|
||||
pokeIndex<LorentzIndex>(Umu,U,nn);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
RealD mass=0.5;
|
||||
RealD M5=1.8;
|
||||
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> HermIndefOp(Ddwf);
|
||||
|
||||
const int nbasis = 8;
|
||||
std::vector<LatticeFermion> subspace(nbasis,FGrid);
|
||||
LatticeFermion noise(FGrid);
|
||||
LatticeFermion ms(FGrid);
|
||||
for(int b=0;b<nbasis;b++){
|
||||
Gamma g5(Gamma::Gamma5);
|
||||
gaussian(RNG5,noise);
|
||||
RealD scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
|
||||
HermIndefOp.HermOp(noise,ms); std::cout << "Noise "<<b<<" Ms "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
|
||||
|
||||
// filter(HermIndefOp,noise,subspace[b]);
|
||||
// inverse iteration
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermDefOp(Ddwf);
|
||||
ConjugateGradient<LatticeFermion> CG(1.0e-5,10000);
|
||||
|
||||
for(int i=0;i<4;i++){
|
||||
|
||||
CG(HermDefOp,noise,subspace[b]);
|
||||
noise = subspace[b];
|
||||
|
||||
scale = pow(norm2(noise),-0.5);
|
||||
noise=noise*scale;
|
||||
HermDefOp.HermOp(noise,ms); std::cout << "filt "<<b<<" <u|H|u> "<<norm2(ms)<< " "<< norm2(noise)<<std::endl;
|
||||
}
|
||||
|
||||
subspace[b] = noise;
|
||||
HermIndefOp.HermOp(subspace[b],ms); std::cout << "Filtered "<<b<<" Ms "<<norm2(ms)<< " "<<norm2(subspace[b]) <<std::endl;
|
||||
}
|
||||
std::cout << "Computed randoms"<< std::endl;
|
||||
|
||||
typedef CoarsenedMatrix<vSpinColourVector,vTComplex,nbasis> LittleDiracOperator;
|
||||
typedef LittleDiracOperator::CoarseVector CoarseVector;
|
||||
|
||||
LittleDiracOperator LittleDiracOp(*Coarse5d);
|
||||
LittleDiracOp.CoarsenOperator(FGrid,HermIndefOp,subspace);
|
||||
|
||||
CoarseVector c_src (Coarse5d);
|
||||
CoarseVector c_res (Coarse5d);
|
||||
gaussian(CRNG,c_src);
|
||||
c_res=zero;
|
||||
|
||||
std::cout << "Solving MCR on coarse space "<< std::endl;
|
||||
HermitianLinearOperator<LittleDiracOperator,CoarseVector> HermIndefLdop(LittleDiracOp);
|
||||
ConjugateResidual<CoarseVector> MCR(1.0e-8,10000);
|
||||
MCR(HermIndefLdop,c_src,c_res);
|
||||
|
||||
std::cout << "Done "<< std::endl;
|
||||
Grid_finalize();
|
||||
}
|
||||
87
tests/Test_cf_coarsen_support.cc
Normal file
87
tests/Test_cf_coarsen_support.cc
Normal file
@ -0,0 +1,87 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=9;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeFermion ref(FGrid); ref=zero;
|
||||
LatticeFermion tmp(FGrid);
|
||||
LatticeFermion err(FGrid);
|
||||
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
|
||||
}
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.8;
|
||||
|
||||
{
|
||||
OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermIndefOp(Dcf);
|
||||
|
||||
HermIndefOp.Op(src,ref);
|
||||
HermIndefOp.OpDiag(src,result);
|
||||
|
||||
for(int d=0;d<4;d++){
|
||||
HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
}
|
||||
err = result-ref;
|
||||
std::cout<<"Error "<<norm2(err)<<std::endl;
|
||||
}
|
||||
|
||||
{
|
||||
OverlapWilsonPartialFractionTanhFermion Dpf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
HermitianLinearOperator<OverlapWilsonPartialFractionTanhFermion,LatticeFermion> HermIndefOp(Dpf);
|
||||
|
||||
HermIndefOp.Op(src,ref);
|
||||
HermIndefOp.OpDiag(src,result);
|
||||
|
||||
for(int d=0;d<4;d++){
|
||||
HermIndefOp.OpDir(src,tmp,d,+1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
HermIndefOp.OpDir(src,tmp,d,-1); result=result+tmp;
|
||||
std::cout<<"dir "<<d<<" tmp "<<norm2(tmp)<<std::endl;
|
||||
}
|
||||
|
||||
err = result-ref;
|
||||
std::cout<<"Error "<<norm2(err)<<std::endl;
|
||||
}
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
58
tests/Test_cf_cr_unprec.cc
Normal file
58
tests/Test_cf_cr_unprec.cc
Normal file
@ -0,0 +1,58 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
const int Ls=9;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
|
||||
}
|
||||
|
||||
RealD mass=0.1;
|
||||
RealD M5=1.8;
|
||||
|
||||
OverlapWilsonContFracTanhFermion Dcf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5,1.0);
|
||||
|
||||
ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
|
||||
|
||||
MdagMLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermPosDefOp(Dcf);
|
||||
MCR(HermPosDefOp,src,result);
|
||||
|
||||
HermitianLinearOperator<OverlapWilsonContFracTanhFermion,LatticeFermion> HermIndefOp(Dcf);
|
||||
MCR(HermIndefOp,src,result);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -34,7 +34,6 @@ int main (int argc, char ** argv)
|
||||
}
|
||||
|
||||
|
||||
|
||||
TComplex cm;
|
||||
|
||||
for(int dir=0;dir<4;dir++){
|
||||
|
||||
63
tests/Test_dwf_cr_unprec.cc
Normal file
63
tests/Test_dwf_cr_unprec.cc
Normal file
@ -0,0 +1,63 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
|
||||
|
||||
const int Ls=8;
|
||||
|
||||
GridCartesian * UGrid = SpaceTimeGrid::makeFourDimGrid(GridDefaultLatt(), GridDefaultSimd(Nd,vComplexF::Nsimd()),GridDefaultMpi());
|
||||
GridRedBlackCartesian * UrbGrid = SpaceTimeGrid::makeFourDimRedBlackGrid(UGrid);
|
||||
GridCartesian * FGrid = SpaceTimeGrid::makeFiveDimGrid(Ls,UGrid);
|
||||
GridRedBlackCartesian * FrbGrid = SpaceTimeGrid::makeFiveDimRedBlackGrid(Ls,UGrid);
|
||||
|
||||
|
||||
std::vector<int> seeds4({1,2,3,4});
|
||||
std::vector<int> seeds5({5,6,7,8});
|
||||
GridParallelRNG RNG5(FGrid); RNG5.SeedFixedIntegers(seeds5);
|
||||
GridParallelRNG RNG4(UGrid); RNG4.SeedFixedIntegers(seeds4);
|
||||
|
||||
LatticeFermion src(FGrid); random(RNG5,src);
|
||||
LatticeFermion result(FGrid); result=zero;
|
||||
LatticeGaugeField Umu(UGrid); random(RNG4,Umu);
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,UGrid);
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
|
||||
}
|
||||
|
||||
ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
|
||||
|
||||
RealD mass=0.5;
|
||||
RealD M5=1.8;
|
||||
DomainWallFermion Ddwf(Umu,*FGrid,*FrbGrid,*UGrid,*UrbGrid,mass,M5);
|
||||
|
||||
MdagMLinearOperator<DomainWallFermion,LatticeFermion> HermOp(Ddwf);
|
||||
MCR(HermOp,src,result);
|
||||
|
||||
Gamma5R5HermitianLinearOperator<DomainWallFermion,LatticeFermion> g5HermOp(Ddwf);
|
||||
MCR(g5HermOp,src,result);
|
||||
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
@ -223,6 +223,9 @@ int main (int argc, char ** argv)
|
||||
//TComplex tracecm= trace(cm);
|
||||
//std::cout << cm << " "<< tracecm << std::endl;
|
||||
|
||||
cm = ProjectOnGroup(cm);
|
||||
|
||||
cm = Exponentiate(cm, 1.0, 10);
|
||||
|
||||
// Foo = Foo+scalar; // LatticeColourMatrix+Scalar
|
||||
// Foo = Foo*scalar; // LatticeColourMatrix*Scalar
|
||||
|
||||
@ -20,12 +20,20 @@ public:
|
||||
sqrtscale = sqrtscale * adj(sqrtscale);// force real pos def
|
||||
scale = sqrtscale * sqrtscale;
|
||||
}
|
||||
// Support for coarsening to a multigrid
|
||||
void OpDiag (const Field &in, Field &out) {};
|
||||
void OpDir (const Field &in, Field &out,int dir,int disp){};
|
||||
|
||||
void Op (const Field &in, Field &out){
|
||||
out = scale * in;
|
||||
}
|
||||
void AdjOp (const Field &in, Field &out){
|
||||
out = scale * in;
|
||||
}
|
||||
void HermOp(const Field &in, Field &out){
|
||||
double n1, n2;
|
||||
HermOpAndNorm(in,out,n1,n2);
|
||||
}
|
||||
void HermOpAndNorm(const Field &in, Field &out,double &n1,double &n2){
|
||||
ComplexD dot;
|
||||
|
||||
|
||||
@ -25,7 +25,6 @@ int main (int argc, char ** argv)
|
||||
std::vector<LatticeColourMatrix> U(4,&Fine);
|
||||
|
||||
NerscField header;
|
||||
|
||||
std::string file("./ckpoint_lat.4000");
|
||||
readNerscConfiguration(Umu,header,file);
|
||||
|
||||
@ -83,7 +82,7 @@ int main (int argc, char ** argv)
|
||||
Complex l = TensorRemove(Tl);
|
||||
std::cout << "calculated link trace " <<l*LinkTraceScale<<std::endl;
|
||||
|
||||
sumBlocks(cPlaq,Plaq);
|
||||
blockSum(cPlaq,Plaq);
|
||||
TComplex TcP = sum(cPlaq);
|
||||
Complex ll= TensorRemove(TcP);
|
||||
std::cout << "coarsened plaquettes sum to " <<ll*PlaqScale<<std::endl;
|
||||
|
||||
@ -39,8 +39,8 @@ int main (int argc, char ** argv)
|
||||
InvSqrt.gnuplot(gnuplot);
|
||||
|
||||
double x=0.6789;
|
||||
double sx=sqrt(x);
|
||||
double ssx=sqrt(sx);
|
||||
double sx=std::sqrt(x);
|
||||
double ssx=std::sqrt(sx);
|
||||
double isx=1.0/sx;
|
||||
double issx=1.0/ssx;
|
||||
|
||||
|
||||
59
tests/Test_wilson_cr_unprec.cc
Normal file
59
tests/Test_wilson_cr_unprec.cc
Normal file
@ -0,0 +1,59 @@
|
||||
#include <Grid.h>
|
||||
|
||||
using namespace std;
|
||||
using namespace Grid;
|
||||
using namespace Grid::QCD;
|
||||
|
||||
template<class d>
|
||||
struct scal {
|
||||
d internal;
|
||||
};
|
||||
|
||||
Gamma::GammaMatrix Gmu [] = {
|
||||
Gamma::GammaX,
|
||||
Gamma::GammaY,
|
||||
Gamma::GammaZ,
|
||||
Gamma::GammaT
|
||||
};
|
||||
|
||||
int main (int argc, char ** argv)
|
||||
{
|
||||
Grid_init(&argc,&argv);
|
||||
|
||||
|
||||
std::vector<int> latt_size = GridDefaultLatt();
|
||||
std::vector<int> simd_layout = GridDefaultSimd(Nd,vComplexF::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);
|
||||
|
||||
LatticeFermion src(&Grid); random(pRNG,src);
|
||||
RealD nrm = norm2(src);
|
||||
LatticeFermion result(&Grid); result=zero;
|
||||
LatticeGaugeField Umu(&Grid); random(pRNG,Umu);
|
||||
|
||||
std::vector<LatticeColourMatrix> U(4,&Grid);
|
||||
|
||||
double volume=1;
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
volume=volume*latt_size[mu];
|
||||
}
|
||||
|
||||
for(int mu=0;mu<Nd;mu++){
|
||||
U[mu] = peekIndex<LorentzIndex>(Umu,mu);
|
||||
}
|
||||
|
||||
RealD mass=0.5;
|
||||
WilsonFermion Dw(Umu,Grid,RBGrid,mass);
|
||||
|
||||
MdagMLinearOperator<WilsonFermion,LatticeFermion> HermOp(Dw);
|
||||
|
||||
ConjugateResidual<LatticeFermion> MCR(1.0e-8,10000);
|
||||
|
||||
MCR(HermOp,src,result);
|
||||
|
||||
Grid_finalize();
|
||||
}
|
||||
Reference in New Issue
Block a user